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INTRODUCTION AND SUMMARY 

In the 1950's, Press and others developed expressions for 
the mean exceedance rate N(y) of an arbitrary aircraft response 
variable through a generic response level y. In deriving these 
expressions, they assumed that the aircraft response Is a locally 
stationary, locally Gaussian random process; the results were 
based on Rice's famous formula for the mean rate of threshold 
crossings of a stationary Gaussian process. Modeling the tur- 
bulence as a locally stationary, locally Gaussian random 
process — generally with a Dryden spectrum — permitted the 
turbulence to be characterized by Its Integral scale and the 
probability density function of Its standard deviation. The 
fact that the standard deviation of the turbulence was treated 
as a slowly fluctuating random variable permitted the mathe- 
matical representation of the turbulence to model the patchllke 
character of real turbulence. 

Measurements recently obtained In a project being 
carried out at the NASA Langley Research Center have demon- 
strated the existence of a low- frequency (large wavenumber) 
component In many turbulence recordings, where this "slow" 
turbulence component Wg(t) appears to fluctuate Independently 
of the patch-llke character associated with the turbulence 
model used by Press and more recent investigators. The addition 
of this large wavenumber component suggests that turbulence 
velocity records w(t) be modeled by a three-component random 
process! 

w(t) = Wg,(t) + w^(t) = Wg(t) ,+ q^(t) z(t), (1.1, 2.3) 

where E{z} = 0, E{z^} = 1 (1.2, 2.5) 

and where the standard deviation of Of(t) of the "fast" turbulence 
component wf(t) satisfies cf(t) ^0. In the work reported here- 
in, we have assumed that the random processes {w (t)}, {cf(t)}, 

s 


"^Most of the equations in this Introductory section have two 
numbers . The first number designates the order of appearance 
of the equation in the present section; the second designates 
the number associated with the same equation, as It appears 
later In the report where the material is treated In detail. 


and {z(t)} ai*e all stationary and mutually Independent. Also, 
we have assumed that z(t) Is a Gaussian process and. In some 
places, that Wg(t) also Is Gaussian. Further discussion of 
this model Is provided In Sec. 2 of the report. 

One of the central tasks of the present work has been to 
determine the conditions that must be satisfied for validity of 
the locally stationary, locally Gaussian response approximation. 
To accomplish this task, we have used the concept of the tur- 
bulence process {w(t)} conditioned on the behavior of the process 
af(t). This conditioning operation Is equivalent to dealing 
with the stochastic behavior of {w(t)} while assuming that the 
function af(t) Is completely specified. Since the processes 
{wg(t)}, {of(t)}, and {z(t)} are assumed to be mutually Inde- 
pendent, this conditioning operation presents no conceptual 
difficulties. Thus, we are able to form expressions for the 
conditional Instantaneous autocorrelation function and Its 
Fourier transform, the conditional instantaneous spectrum of 
the turbulence process {w(t)}, given that the function af(t) 

Is specified. These expressions are derived in Sec. 3-1- 

Although the process {w(t)} is stationary, the process 
{w(t)} conditioned on af(t) is. In general, nonstationary. 
However, in an earlier study that dealt with the effects of 
nonstationary behavior on the spectra of atmospheric turbulence, 
a series expansion of the instantaneous spectrum was developed 
for studying the effects of the time variations of Of(t) on the 
Instantaneous spectrum. Some results of this earlier study , 
relevant to the present work, are reviewed In Sec. 3*2. The 
first term in this series expansion, when applied to the "fast" 
component Wf(t) In Eq. (1.1)' and conditioned on the function 
af(t). Is the usual locally stationary spectrum approximation 


(f,t|a^) s a^(t) 


(1.3, 3.11) 


where 'J’zCf) is the power spectrum of the stationary process 
{z(t)} of Eq . (1.1). Thus, Investigation of the correction 
terms provided by the series expansion have enabled us to formu- 
late conditions for validity of the locally stationary response 
approximation. 

First, it was necessary to derive a series expansion for 
the Instantaneous spectrum of the response process, also condi- 
tioned on the behavior of ap(t). To obtain this expansion, we 
have used the Input-response relationship for the Instantaneous 
spectrum derived In our earlier report. The conditioned Instan- 
taneous spectrum of the aircraft response process {y(t)> is 
given by Eq . (3.29) of the present report. When the series 
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expansion for I cff ) is combined with this Input-response 

relationship, we obtain the series expansion for the conditioned 
Instantaneous response spectrum given by Eq. (3* 37) or ( 3 . 38 ) 
and derived In Sec. 3-^. The leading terms In this response 
representation are 

$y(f,tlaf) ~ (f)+aj(t)$^(f)] |H(f)|2, (1.4, 3.40) 

s 

where (f) Is the power spectrum of the "slow" component Wg(t) 

In the turbulence model of Eq. (1.1) and H(f) is the aircraft 
complex frequency-response function. Equation (1.4) is the 
obvious locally stationary response approximation that could 
have been written directly from Eq. (1.1). 

Examination of the appropriate correction terms In 
Eq. ( 3 * 38 ) to the locally stationary response approximation 
given by Eq. (1.4) has enabled us to write conditions for 
validity of the locally stationary response approximation of 
Eq. (1.4) in Sec. 3.5. Three conditions, Eqs . (3.41), (3-43)j 
and (3.46) are given. Equation (3.41) expressed the local 
stationarlty requirement for the turbulence w(t) of Eq . (1.1), 
whereas Eqs. (3-43) and (3-46) express the local stationarlty 
requirements for the aircraft response y(t), assuming that the 
requirement of Eq. (3.41) Is satisfied. 

The local stationarlty requirements of Eqs. (3-41), (3.43)* 
and (3-46) are expressed^ln terms of (derivatives of the loga- 
rithm of) the function af(t), which Is still assumed to be a 
specified function at this juncture. Before discussing these 
requirements and formulating them In terms of stochastic 
metrics of Of(t), we derive expressions for the aircraft- 
response exceedance rate N+(y) and the first-order probability 
density function p(y). These expressions are derived In Sec. 4, 
where we assume validity of the locally stationary response 
approximation given by Eq . (1.4). 

In Sec. 4.1, It Is shown that If the process {z(t)> In 
the model of Eq . (1.1) Is Gaussian then the process- {w.f'(t ) 1 Cf } , 
conditioned on the process Cf(t), also is Gaussian — this 
result being Independent of locally stationary requirements. 
However, If probability density functions of w^(t) are generated 
by time averages, then local stationarlty Is required for the 
"fast" process (wf(t)} to be considered locally Gaussian. 

In Sec. 4.2, we derive the expression for response exceed- 
ance rates (with positive slopes) given by 
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0 


(1.5, 4.8) 


N^(y) 


N^(y 


op 


p(apda^ , 


where p(a|>) is th^ probability density of the square of the 
process af>(t) — af(t), therefore, being the local variance of 
the fast turbulence component Wf(t) — and where N+(y|af) is 
the expected local up-crossing rate through the threshold y 
of the response process, given that the local value of af(t) 
is specified and assuming the locally stationary response 
approximation of Eq. (1.4) to be justified. The expression 
for N 4 .(y|af) is given by Eq. (4.22). To evaluate N+(y|ap we 
require spectra $g;(^) (^) 4he turbulence components 

s . 

z(t) and Wg(t), as well as the magnitude of .the aircraft 
frequency-response function.. 

For cases where va:^latlons In Cf(t) are small relative to 
the expected value of ap(t), a useful series expansion for 
N+(y) has been derived in Sec. 4.3. This series expansion is 
of the form 


oo . 

N+(y) = I ’ 

k=0 ■ °f 


where we have defined 



(1.6, 4.28) 


(1.7, 4.26) 


where 


= E{op 


(1.8, 4.24) 


is the expected value of and where 




(apop 


k 


P(ap 


d(op 


k = 1,2, 


(1.9, 4.29) 


are the central moments of apt). Of particular Interest is the 
two-term approximation of Eq. (1.6) given by 
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N^(y) ~ N^(y|ap 


+ 


( 1 . 10 , 


^. 31 ) 




N|^^(ylop 


= N^(y|ap [l+y^|^Q^^^y|op], ( 1 . 11 , 4.36) 

where N^^^(y|ap has been written in the second line as 

= 2N^(y|^) Q^^^(yl^)- d-i2, 4.33) 

From Eq. (1.10), It Is evident that the first term N^(y]ap 
on the right-hand side Is the exceedance rate one obtains by 
assuming that y(t) Is a stationary Gaussian process.* Thus, 
the second term In Eq. (1.10) is a correction term that modifies 

the Gaussian approximation N+(y(cr|.) to account for relatively 
small fluctuations In o?.(t). When written In the form of 

•T ( 2 ) - .y 

Eq. (1.11), the general form of the correction term (yjop 
Is given by Eq. (4.34) In the text. 


For cases where the response of the aircraft to the slow 
component Wg(t) is negligible in comparison with the response 
to the component Cf(t) z(t), the form of N+(y) given by Eq. 
(1.11) is particularly instructive. For these cases, when we 


•Equations (1.6) and (1.10) provide motivation for expressing 
N^(y) as an integral over the probability density p(cr^) of the 
local variance rather than as an integral over the proba- 
bility density or the local standard deviation Of, as was 
done by Press and others. The first term in the series 
expansion of N^(y) given by Eq. (1.6) is, according to 

Eq. (1.10), N (y|^). However, N (y|cl) is simply the 
exceedance rate obtained by assuming tnat a^-Ct) Is a constant. 
In cases where 0 f(t) varies somewhat with time and we estimate 
N+(y) by assuming that the response is a stationary Gaussian 
process, we obtain for our estimate of N+(y) the first term 

N^(y|op 'the right-hand sides of Eqs. (1.6) and (1.10) as 

is shown in Sec. 4.3 [see Eq. (4.32)]. However, since 

E{ap 7 ^ {E[a^]}^, the quantity N^(y|o^) is different from the 

exceedance rate that would be obtained by evaluating the 
expression for the exceedance rate for stationary Gaussian 
processes using for the mean of the probability density of 

rather than the square root of 
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( 2 ) 

take the logarithm of Eq..(l.ll) and assume that y^ 2 '^is small 

enough so that we may use £,n(l+x) = x, we obtain In Sec. 4.3 
the simple relationship 


£n 


N^(y) 


N 


( 0 ) 


. _ JLl. 


1 




(ap^ 




(1.13, 4.40) 


The first term (-y^/2ap on the right-hand side of Eq. 

(1.13) Is the familiar result for stationary, Gaussian processes; 
l.e., the logarithm of the normalized exceedance rate 
N+(y)/N+(0) Is linear when plotted as a function of the square 
of the response level, and is the variance of the response 


y 


process. The coefficient 






governs the strength of 


the correction to the Gaussian approximation given by the first 
term. The quantity , V( ap ^ Is the square of the coefficient 

U ^ J- 

of variation of the time-varying variance apt) In the model of 

( 2 ) 

Eq. (1.1). Except for the value of y 2 V(al)^, which is always 

I 

nonnegative, the functional form of the correction term In the 
right-hand side of Eq. (1.13) is fixed. For y = 0, the correc- 
tion Is zero; hence, the Gaussian approximation given by the 
first term yields the correct value. For 0 < y < (2a ), the 
correction term Is negative; hence, in this Interval, ^the 
Gaussian approximation overestimates the threshold-crossing 
rate. For y > ( 20 y), the correction term Is positive; hence, 
for large values or y the Gaussian approximation underestimates 
the threshold-crossing rate. These theoretical predictions are 
consistent with well known experimental observations, when 
considered as a function of y^ (or y^/a| ), Eq. (I. 13 ) predicts 
(a parabolic) concave exceedance plot. 


A treatment of the probability density function p(y) of 
the aircraft response Is provided In Sec. 4.4, where we have 
used the same general approach that was used for exceedance 
rates. For example, it is shown there that we may express 
P(y) by 


p(y) 


p(y I app(ap da^ . 
0 


(1.14, 4.42) 
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A series expansion for p(y) that Is developed [Eq. (U.US)] Is 
analogous to that of Eq. (1.6). The two-term approximation to 
this series expansion that Is analogous to Eq. (1.1) Is 


p(y) = p(y|ap 
where U^^^(ylap Is 

U^^^yl'^) = I 


[l+y^f^U^2\y|ap] 

Of I 


defined by 



(1.15, ^.51) 


(1.16, 4.50) 


where we have used the notation of Eq. (4.48). The quantity 
of Is the aircraft mean-square response to turbulence component 
z 

z(t) given by Eq. (4.l8), and Oy Is the overall mean-square 
response given by Eqs. (4.l6) and (4.48). It Is shown In 
In Sec. 4.4 that the correction term given In Eq . (1.15) to 

the Gaussian approximation p(y|ap of p(y) Is exactly the same 
as the first correction term provided by the Gram-Charller 
expansion. However, the derivation of Eq. (1.15) Is based on 
what would appear to be a completely different line of rea- 
soning. 


For cases where the aircraft response to the "slow" tur- 
bulence component w„(t) Is negligible, the correction term to 
the Gaussian approximation In Eq. (1.11) Is given by 


Q 


( 2 ) 




(1.17, 4.39) 


Comparison of Eqs. (I.l6) and (1.17) Indicates that whenever 
the fluctuation In apt) Is not negligible In comparison with 

its mean apt), the first-order probability density p(y) and 
the exceedance rate N+(y) have different functional dependencies 

( 2 ) 

on y. However, when y 2 = 0, both p(y) and N+(y) have the 

Of 

shape of Gaussian probability density functions. 

The 'locally stationary response requirements. Initially , 
treated In Sec. 3, are discussed and related to stochastic 
metrics of the process a^(t) In Sec. 5. In particular, it is 
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shown In Sec. 5-1 that the locally stationary response require- 
ments may be Interpreted as requiring (1) that the relative 
changes In af(t) must occur slowly when measured on the time 
scale t' = L 2 /V, where Is the Integral scale of the component 
z(t) In the model of Eq. (1.1) and V Is the aircraft speed 
[see Eq. (5-1)]; (2) that the relative changes In af(t) must be 
small In comparison with unity when measured over time Intervals 
of the order of the group delay of the aircraft Impulse response 
function [se| Eqs . (5.^) and (5-6)]; and (3) that the relative 
changes In af(t) must occur slowly relative to time intervals 
comparable to the nominal duration of h^(t), where h(t) Is the 
aircraft unit-impulse response function [see Eqs. (3.^6) and 
(5.7)]. 

These three spectrum conditions are expressed In terms of 
(somewhat less stringent) mean-square response conditions by 
Eqs. (5.8)» (5.9)j and (5.10) in Sec. 5.2. In Sec. 5.3, the 
requirements on the behavior of a^(t) are expressed In terms 
of derivatives of the autocorrelation function of the logarithm 
of af(t),^where the dependence on the logarithm of af(t), rather 
than on af«(t) Itself, has been dictated by the requirements 
themselves. [It Is the fvaational fluctuation of Of(t) that 
is Important, rather than the absolute fluctuation of a|-(t)]. 

The final forms of the three local statlonarity requirements 
are* 


[r|^^0)]'^ < 3.27t2 


$^(f) |H(f) l^df 


— 00 


(f ) |H(f ) I ==df 

— 00 


[-r;' (0)]'^ 



j $^(f)mj^°^(f)df 

— 00 

»oo 
— 00 


( 1 . 18 , 5 . 17 ) 


( 1 . 19 , 5 . 18 ) 


*0ne of the other forms of these requirements given In the main 
text In Sec. 3 or 5 may be easier to apply In practice. Equa- 
tions ( 1 . 18 ), ( 1 . 19 ), and (1.20) would seem to be the least 
restrictive set of requirements. 
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and 


- i R^'(O) (l + 3 { 2 


1 »l5 


[r; ’( o)]^( 


-oo . 

J (f)df 


$^(f)m^^^(f)df 


( 1 . 20 , 5 . 22 ) 


,( 4 ) 


where Ry'(O) and '^(O) are the second and fourth devlatlves 
of the autocorrelal^ion function of 


v(t) = ilna^Ct), (1.21, 5.11) 

where these derivatives are to be evaluated at the origin, $ (f) 
is the power spectrum of the component z(t) In the model of ^ 

Eq. (1.1), $^^^(f) Is the second derivative of H(f) 

is the aircraft complex frequency-response function, and for 

n = 0, 1, and 2, m;^'^^(f) is the power-moment spectrum of the 
aircraft unit-impulse response function h(t) defined in earlier 
work by this writer in a completely different context: 



(f) 


4 


«00 

t^$^(f ,t )dt , 

— oo 


( 1 . 22 , 3 . 36 ) 


where $^(f,t) is the Instantaneous spectrum of h(t) defined by 
Eqs. ( 3 . 25 ) and (3.26). Writing the complex frequency-response 
function H(f) as 


H(f) 


H(f) le 


10h(f) 


(1.23, 5.3) 


we express m^^^(f), m^^^(f), and m^^^(f) in terms of the magni- 
tude and phase of H(f) by 


m^°^f) 


= |H(f)| 


m^^^f) = - ^ |H(f) 


d0, (f) 

2 h 
df 


(1.24, 3 . 39 ) 
(1.25) 
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1 


( 1 . 26 ) 



(f) 


47T‘ 


|H(f) 



1 d^&n|H(f ) I ) 

2 df^ ] ’ 


where the time origin of h(t) must be chosen to satisfy Eq. 

(3.22) and where Eqs. (1.25) and (1.26) are a consequence of 
Eqs. 03.39), (5.^), and (5.7). Equations (l.l8), (1.19), and 
(1.20) would appear to be the least restrictive requirements 
that must be satisfied for confident engineering usage of 
the locally stationary response approximation of Eq. (1.4), 
which is the basic approximation used in deriving the above 
described expressions for the exceedance rates and probability 
density functions of an aircraft-response variable y(t). 

Evaluation of the left-hand side of the requirements of 
Eqs. (1.18), (1.19), and (1.20) requires a capability to evaluate 
from measured turbulence records the autocorrelation function 
Ry(T) of ZnOf(t), as indicated by Eq. (1.21). Since only the 
derivatives of Ry(T) are required, Ry(t) need be determined 
only to within an additive constant. In Sec. 5.4, it is shown 
that Rv(t) may be expressed as 


R^(t) = {E[iln a|(t)]}2 + Cov[£n w^(t), Hn w^(t+r)] 

- 2 arcsln^ [R (t)/R^ (0)], (1.27, 5-36) 

^h ^h 


where {E[£n Cf(t)]}^ is a constant, Cov[***] is the covariance 
of An w^(t) and An w^(t+x), which can be evaluated directly 
from a high-pass filtered version w^(t) of the turbulence 
record w(t) as described in Sec. 5.4, and R (t) Is the (in- 

^h 

verse) Fourier transform of the high-pass filtered version of 
$g(f) defined by Eq. (5.28). All of the above quantities are 
amenable to numerical calculation from turbulence velocity 
records. To obtain Eq. (1.27), we had to derive an expression 
for the autocorrelation function of the natural logarithm of 
the square of a stationary Gaussian process from the autocorre- 
lation function of the process Itself. This result is given 
(for processes with zero mean and unit variance) by Eq. (5.31). 

The above work suggests that, in addition to the auto- 
correlation function Rv(t) of An a|(t), an adequate characteriza- 
tion of turbulence records for aircraft response predictions 

requires the spectra $ (f) and $ (f) of the components w (t) 

Wg z b 

and z(t) in the turbulence model of Eq. (1.1) and the probability 
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density function p(af) of the square of the component a^Ct) In 
Eq. (1.1). All of these quantities are required for computation 
of response exceedance rates - e.g., see Eqs . (1.^) and (1.5). 

In addition, it would be desirable to test the assumption that 
{Wg(t)> is a (stationary) Gaussian process. The "first-order 
test" of this assumption is the probability density function 
of the component Wg(t). Especially relevant metrics of the 
probability density functions of oMt) and Wg(t) are the low- 
order central moments — e.g., see Eq. (1.13). Furthermore, 
it appears likely that the fourth-order moment, of the response 
process y(t) can be computed if the power spectrum of a fit) is 
available, even in cases where neither of the locally stationary 
response conditions of Eqs. (1.19) and (1.20) are satisfied, 
but where the requirement of Eq. (I.l8) is satisfied.* Thus, in 
addition to the autocorrelation function of Jin a|>(t), the power 
spectra of z(t), Wg(t) , and crf(t) and the moments and prob- 
ability density functions of Of(t) and Wg(t) are useful turbu- 
lence characterizations for the prediction of aircraft-response 
statistics. Methods for computation of these turbulence metrics 
are developed in Sec. 6. Specifically, methods for estimating 
the spectra of Ws(t) and z(t) ar| described in Sec. 6.1, methods 
for estimating the spectrum of Cf(t) are described in Sec. 6.2, 
and methods foy estimating the moments and probability density 
functions of Of(t) and Wg(t) are described, respectively, in 
Secs. 6.3 and 6 . . 

As is the case with most newly developed research results, 
the methods and conclusions reported herein represent a (hope- 
fully error free) "first cut" at the problem of adequately 
taking into account — in a not too complicated fashion — the 
nonGausslan behavior of real turbulence records for the purpose 
of predicting aircraft-response statistics. Recommended future 
work would include "fine tuning" and improvements to methods 
and conclusions reported herein. 

Aaknowledgements . Comments by Dr. Kenneth Sidwell have 
been helpful in several stages of this work as detailed in the 
footnotes. The typing of the report was carried out by Mrs. 

Linda Nelson, and the figures were prepared by Ms. Marla Maldter. 
The support of the work by NASA Langley Research Center is grate- 
fully acknowledged. 


*This result will be of use in the case of very high-speed 
aircraft . 
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NONGAUSSIAN TURBULENCE MODEL 


Atmospheric turbulence velocity measurements are usually 
made In sets of three records: the vertical, lateral, and 

longitudinal time histories. When Taylor's hypothesis — l.e., 

X =. Vt - Is employed, there result vertical, lateral, and longl 
tudlnal records that are regarded as functions of a spatial 
variable x. The "standard" turbulence model Is to assume that 
each of these records Is a sample function drawn from a station 
ary (or homogeneous) Gaussian random process. Furthermore, It 
Is generally assumed that the power spectra of the vertical and 
lateral records should be well described by the von Karman 
transverse spectrum and that the longitudinal records should 
be well described by the von Karman longitudinal spectrum. 

Each of these two von Karman spectral forms Is described by two 
parameters, the rms turbulence velocity a and the Integral 
scale L. Since the power spectral density provides a complete 
probabilistic description of a Gaussian random process, these 
von Karman spectral forms provide complete statistical descrip- 
tions of the three components of measured turbulence veloci- 
ties — If the stationary Gaussian model with von Karman spectra 
is a valid model of turbulence velocity records. 

However, results predicted by the above model are not 
always In agreement with turbulence measurements. Differences 
between the model and observed records are manifested in at 
least two ways: (1) the probabilities of large excursions are 

underestimated by the model, and (2) the low-frequency content 
of turbulence velocity spectra Is often underestimated by the 
model. In addition, the "knee" of measured spectra Is often 
less sharp than the knee of the appropriate von Karman spectral 
form. 


The fact that the probability of large excursions Is often 
underestimated by the standard turbulence model Is directly 
attributable to the fact that turbulence time histories often 
appear to have a time-varying envelope or patchy character. 

One may regard such time histories to be either nonstationary 
but Gaussian or stationary but nonGaussian. * For a given 
standard deviation (rms value) of turbulence velocity, a time 
variation In the envelope generally has the effect of yielding 
more large excursions than the Gaussian probability density 
predicts. This behavior accounts for the "concave shape" of 
turbulence exceedance plots. 


*The distinction between nonstationary Gaussian behavior and 
stationary nonGaussian behavior Is, In general, impossible to 
make with a single time history. It Is often more convenient 
to regard such records to be stationary but nonGaussian (or 
locally Gaussian). 
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The above comments suggest that a given turbulence record, 
say w(t), may be modeled by , 

w(t)=a(t)z(t), a(t)>0 (2.1) 

E{z} = 0, E{z*} = :i, (2.2) 

where a(t) may be regarded as a time-varying standard deviation 
(or envelope), and where {z(t)} is a stationary process, with 
zero mean value and unit standard deviation, which may be taken 
to be Gaussian. Time variations of the function a(t) account 
for the patch-like character of atmospheric turbulence veloci- 
ties observed in many records. Insofar, as the visual appearance 
of the velocity records is concerned, the function a(t) may be 
regarded as either a stochastic dr a deterministic function. 
Appropriate choice of the amplitude distribution of o(t) will 
allow the probability distribution of w(t) to take on a wide 
variety of forms. 

Visual inspection of turbulence records indicates one addi- 
tional feature that cannot be accounted for by the model of 
Eq. (2.1). Some records which exhibit the patch-like character 
modeled by Eq. (2.1) also have superimposed on this patch-like 
structure another very low frequency component, which appears 
to fluctuate Independently of the envelope of the higher- 
frequency components that possess, the patch-like character. 

This apparent independence suggests that a low-frequency term 


be added to Eq. 

(2.1); i.e. , ■ : 


w(t) 

= Wg(t) + Wj,(t) 


where 

= Wg(t) + a^(t) z(t). 

(2.3) 

Wj.(t ) 

= a^(t) z(t), a^(t) ^ 0 

(2.M 

and 



E{z } 

= 0, E{zM = 1. 

(2.5) 


In the following work, we shall assume that {wg(t)>, (af(t)}, 
and {z(t)> is each a stationary random process and that 
(z(t)} is a Gaussian process with. zero mean value and unit 
variance, as indicated by Eq. (2.5). The spectrum of the 
"slow" process {wg(t)> generally occupies a frequency range 
that is low in comparison with that occupied by the "fast" 
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process {wf(t)}; it Is th.e low-frequency contribution of {Wg(t)} 
that causes underprediction by the von Karman model of the 
low-frequency part of the spectrum. Thus, we shall further 
assume that the power spectrum of {z(t)} is described by the 
appropriate von Karman form and that {wgCt)}, {a^Ct)}, and 
{z(t)} are mutually Independent processes. In some of our 
work, we shall also assume that {wg(t)} is a zero mean Gaussian 
process. However, we shall not want to assume that {a(t)} is 
Gaussian, since such an assumption would permit a(t) to go 
negative. Equation (2.3) is the simplest model of atmospheric 
turbulence that possesses the flexibility required to represent 
readily observable features of measured turbulence records. 

At this Juncture, it is appropriate to illustrate the need 
for the model of Eq. (2.3) by examination of some measured tur- 
bulence velocity histories. The vertical record shown in Fig. 1 
(Ref. 1) Illustrates a record that is probably reasonably well 
modeled by a stationary Gaussian process, especially the 
right-hand half of the record — l.e., from 120 sec elapsed 
time to the end. Thus, this record would be reasonably modeled 
by Eq . (2.3) with Of(t) set equal to a constant — in which case 

there is no need for the term w (t). 

s 

The records shown in Fig. 2 (Ref. 1) Illustrate mild patch- 
like behavior. For example, the patches occurring at 150 and 160 
sec elapsed time on the vertical record illustrate distinctly 
nonstationary or nonGaussian behavior. The number of large 
excursions of the records shown in Fig. 2 is substantially 
larger than would occur for stationary Gaussian records with 
the same standard deviations and spectra. 

Each of the records shown in Fn.g. 2 also exhibits an addi- 
tive low-frequency component that appears to fluctuate indepen- 
dently of the occurrence of the patches. For example, for the 
5-sec Interval between 183 and l88 sec on the vertical record, 
high-frequency fluctuations are absent; however, there remains 
in that Interval a fluctuating low-frequency component. Similar 
behavior occurs between approximately 96 and 99 sec elapsed 
time on the vertical record shown in Fig. 3 (Ref. 1). A strong 
low-frequency component is present there; however, high-frequency 
fluctuations are absent in that interval. These "gaps" cannot 
be explained in terms of statistical fluctuations of a station- 
ary Gaussian process. 

Each of the records shown in Fig. 4 (Ref. 2) dramatically 
Illustrates the three components of the turbulence model of Eq . 
(2.3). For example, in the region from 9 min 0 sec to 9 min 
45 sec on the vertical record, the term a^(t) in Eq . (2.3) 
is negligible in comparison with the very strong low-frequency 
component Wg(t), which is clearly present over the entire 
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record. At about 9 min sec, af(t) grows and then decays 
back to a small value about one minute later. The behavior 
of the records shown in Fig. aannot be modeled by a single 
stationary Gaussian process. 

A model functionally similar to Eq . (2.3) has been sug- 
gested by Reeves et at. [3,4] for use with flight simulators. 
However, in Reeves' model, {o(t)} and { 2 (t)} are specified as 
(stationary) Gaussian processes, and the processes {v/g(t)} and 
iwf(t)} both have the same spectral form — the Dryden spectrum - 
and both have the same integral scale. The model of Eq. (2.1) 
has been proposed independently by Sidwell [ 5 ] and Mark [fi]. In 
addition, Sidwell in work currently being carried out at 
NASA Langley Research Center, has proposed the addition of the 
low-frequency component v^g(t) indicated by the model of Eq . (2.3). 

The model of Eq. (2.3) could be generalized further — e.g., 
by including a multiplicative modulating term ag(t) in the 
slow component Wg(t). In this case, the slow and fast compon- 
ents Wg(t) and w^(t) would have identical functional forms. 

This generalization, which may have to be made after further 
examination of additional turbulence records, requires only 
modest changes in the methods developed in this report. 
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LOCALLY STATIONARY TURBULENCE AND AIRCRAFT RESPONSE APPROXIMATIONS 


Conditional Instantaneous Spectrum of Turbulence Model 


In the turbulence velocity model described by Eq . (2.3) > the 
three stochastic functions Ws(t), of(t), and z(t) are assumed to 
be mutually independent sample functions drawn from the three pro- 
cesses {ws(t)}, {af(t)l, and {zCt)}. For reasons that will become 
evident later, we shall be concerned here with the conditional 
instantaneous spectrum of the stochastic function defined by Eq . 
(2.3); i.e., 

w(t)=Wg(t)+aj.(t)z(t) , (3.1) 

where this instantaneous spectrum will be conditioned on the 
stochastic function Of(t). The operation of forming the conditional 
expectation required in developing the instantaneous spectrum of 
w(t) is equivalent to treating df(t) as a known (or deterministic) 
function of completely general form, except for the requirement 
that af(t) _> 0.* 

To determine an expression for the conditional Instantaneous 
spectrum of w(t), we first form the conditional instantaneous auto- 
correlation function of w(t). The required instantaneous auto- 
correlation function <()w(t jt I Of) is defined — e.g., Mark [7], Mark 
and Fischer [ — as 

4>^^(T,t lo^) = E{w(t-^) w(t+^) 1 o^(t ) } (3.2a) 

= E{[w^(t-^) + a^(t-^) z(t-^)] 

[Wg(t+^) + aj,(t+^) z(t+^) ] |o^(t) } , (3.2b) 

where the vertical bars followed by Of(t) In the left- and right- 
hand sides of Eq. (3-2) indicate that the expectation operation 
E{...|af(t)} assumes that the (stochastic) function ap(t) is known 
or specified; hence, no averaging is carried out over the ensemble 
of functions {op(t)}. See, for example, pp . 55 and 56 of Lanlng 
and Battln [S] for a brief discussion of the notion of a conditional 
expectation. (Chapters 1 through 5 of this reference provide an 
excellent introductory discussion of probability theory and 
stochastic processes.) Expanding the right-hand side of Eq . (3.2b) 


The concept of the process wf(t) conditioned on the process Of(t) 
has been used by Sidwell [5]. 
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and recognizing that {wg(t)} and {z(t)} are independent stationary 
processes, that E{z} = 0, and that C£.(t) is being treated as a 
knov/n function, we have 

(f)^^(T,t loj,) = E{Wg(t-|-) Wg(t+|-)} 

+ a^(t-^) o^(t+^) E{z(t-^) z(t+^)} 

= 4 >„ (t) + (}>^^{T,t lo^) > ( 3 . 3 ) 


where 4>wg(T) and are the autocorrelation functions of {wgCt)} 

and {z(t)} and where we have defined 



(T,t|a^) = a^(t-^) a^(t+|-) , 


( 3 .^) 


which requires no stochastic average since a^(t) is assumed to be 
specified . 

The (conditional) instantaneous spectrum of w(t) is defined 
[e.g., 5,7] as the Fourier transform with respect to t of 
<|)w(T,t jof) : 




<f)^^(T,t [oj.) dT 


( 3 . 5 ) 


Recognizing that the Fourier transform of a product is the con- 
volution of the transforms, and treating t as a parameter in Eq . 
(3.3), the conditional instantaneous spectrum of w(t) may be 
Immediately obtained from Eq . (3-3) as 


$„(f,tla^) = (f) + 


(S,t 


a^) $^(f-g) dg 


( 3 . 6 ) 


where (f) and 4>z(f) are the power spectra of the stationary 
s 

processes {ws(t)) and (z(t)}, which are obtained by forming the 
Fourier transforms of the appropriate autocorrelation functions — 
i . e . , 
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. i N -i2irfT , 
<f>. (t ) e dT 

— OO 


(3.7) 


— and where the conditional instantaneous spectrum 4>Q^(f,tjof) 

in Eq . (3-6) is obtained by forming the Fourier transform with 
respect to t of x , t | ) . Equation (3-6) is the desired ex- 


pression for the conditional Instantaneous spectrum of the pro- 
cess (w(t)} defined by Eq. (3.1). 


Series Expansion 
Spectrum 


of Conditional Instantaneous 
of Turbulence Model 


Equation (3.6) can be written as 

(f) + 4^ (f.t|o^) , 

s f 


where 




|a^) $^(f-g) dg 


(3-8) 


(3.9) 


is the conditional instantaneous spectrum of the component w„(t) 
defined by Eq . (2.H); i.e.. 


w^(t) = a^(t) z(t) . (3.10) 

If the temporal variations in a^(t) occur slov;ly in comparison 
with those of z(t), we v/ould expect the conditional instantaneous 
spectrum of i^r^(t) to have the form 

<f^^^(f,tla^) aj(t)- 4>^(f)- , (3-11) 

where 'J>z(f) is the pov/er spectral density of the process (z(t)}. 
Equation (3.12) is the Zocallij stationary approximation for the 
oonditional instantaneous spectrum of the "fast" component of 
turbulence w^(t). 

For most measured turbulence records, the fluctuations of 
Of(t) do, in fact, vary slowly in comparison v/ith those of z(t). 



Mark and Fischer [&] derived a series representation of the con- 
ditional Instantaneous spectrum f ,t | af>) that enables one to 

determine the conditions required for the approximation of Eq . 

(3- 11) to be valid. This series representation v/as developed in 
terms of distance x and wavenumber k, rather than in terms of 
time t and frequency f. However, the results of Ref. 6 may be 
applied directly to the variables t and f simply by substituting 
t for X and f for k. Applying this substitution, we find from 
Eq. (4,11) of Ref. 6, with minor changes in notation. 


f (f ,t 

l.T ' ^ 


V/ 


O^) 


N 

= I 

n=0 


an(t) 

nl 


(f ) 

z 


R„„pr.t) 


(3.12) 


v/here, from Eq . (4.9) of Ref. 6, (f) is defined as the nth 

derivative of the power spectral density <I> 2 ;(f) of (z(t)) — i.e.. 


$(n)(^) 4 ^ 
^ df 


n 

- $ (f) 
n z ' 


and where 


$(0) (f ) = $ (f ) 

z z 


(3.13) 


(3.14) 


The coefficients an(t) in Eq . (3-12) may be expressed in terms of 

the derivatives of Of(t) by 


^n^^> = TTtrn ^ 

(-i4-^) k=o 


n 

I 


(- 1 )^ 


(*^) 


d'^a^(t) d^ ^a^(t) 


dt 


k 


dt 


n-k 


(3.15) 


according to Eq . (4.19) of Ref. b. The quantities (j^) are the 
binomial coefficients. From the above expression for aj.j(t), one 
may show that for odd integer values of n. 


an(t) - 0 , n = odd 


(3.16) 


Expressions for the remainder term ^^9 • (3*12) are 

given by Eqs. (4.7) and (4.13) of Ref. 6. We are particularly 
interested in the first two nonvanishing terms a^(t) which, 
according to Eqs. (4.50) and (4.52) of Ref, 6, may be expressed 
as 

a,(t) = qj(t) (3.17) 
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aj (t ) 


1 

8tt" 


(3.l8a) 




(t) 


d^Jlna^Ct ) 
dt^ 


16t' 


olit) 


d^iln[a^(t ) ] 
dt* 


(3.l8b) 


Combining Eqs. (3.12), (3-1^), (3.l6), (3.17), and (3.18), we 
have for the first two nonvanishing terms in our series expansion 
for <I'vyj.( f ,t jof ) 




a^(t ) 




32n' 


d^S-no^(t) 

dt^ 


4>2^^{f) 


+ . 


(3.19) 


The first term in Eq. (3.19) is the locally stationary spectrum 
approximation given by Eq . (3.11). Whenever the second term in 
Eq. (3.19) Is negligible In comparison with the first, the locally 
stationary approximation is valid. This condition may be expressed 
as 


JlnOp(t ) 
dt^ 


<< 32 tt^ 


. VIL_ 
1 


( 3 . 20 ) 


The above condition will be expressed in terms of measurable 
turbulence velocity metrics in a later section of this report. 


By combining Eqs. (3.8) and (3.19), we obtain the series 
approximation to the conditional instantaneous spectrum of both 
turbulence components; i.e.. 


\(f,t 


Of) = 


4>„, (f) + oi(t) 


w 


»,<n - 


32n- 


Znoi(t ) 


dt' 


4,(2)(f) 


+ 


3 


( 3 . 21 ) 


which is the desired result. Notice that the correction term to 
the locally stationary approximation involves the second time 
derivative of £nOf(t). Thus, when £na|>Ct) fluctuates sufficiently 
slowly — i.e., when Eq. (3.20) is satisfied — the locally station- 
ary approximation is valid. 
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Conditional Instantaneous Spectrum of Aircraft Response 


Consider an aircraft with spatial dimensions negligible in 
comparison with the scale of the von Karman component z(t) of 
the turbulence model of Eq . (2.3). We may characterize a desired 
aircraft response variable by the.response h(t) of this variable 
to a "unit Impulse" of turbulence velocity w(t) applied at t = 0. 
However, for reasons that will become evident later, we shall 
want to choose the position of the time origin of h(t) so that 
the time centroid of h^(t) occurs at t = 0; l.e., the time ori- 
gin shall be chosen so that h(t) satisfies 


t h^ (t )dt 

. 

..oo 


h^ (t )dt 

_00 


0 


( 3 . 22 ) 


Equation (3.22) defines a unique position for the time origin of 
h(t); it is located at the center-of-gravity of the "mass dis- 
tribution" h^(t). 


Let us denote the conventionally defined unit-impulse 
response function by h(t)j_ where the time origin of h(t) is 
generally chosen so that h(t) = 0 for t < 0. The time centroid 
of h^ (t ), which, we shall denote by tjy, is defined by 


r~ _ 

t h^ (t )dt 



Leo 

»00 

h^ (t )dt 

~oo 


(3.23) 


Our definition of h(t) that_satisfles Eq. (3.22) may be related to 
the conventionally defined h(t) by 

h(t) = h(t+tp^) , • (3.24) 

which can easily be verified by substitution of Eq . (3.24) into 
Eq . _^3.22), Introducing a change of variable, and then, solving 
for tp^. The result of these operations yields Eq . (3.23). 
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To derive the desired’ Input-response relationships for air- 
craft, we define the instantaneous autocorrelation function of : 
the (deterministic) unit-impulse response function h(t) as 


<J>h('r,t) - h(t-^) h(t+^) , 


(3.25) 


from which we may define the Instantaneous spectrum of h(t) as 

( 3 . 26 ) 




4. / -12TTfT , 

(|)j^(T,t)e dT 


Let {y(t)} denote the generally nonstationary response process of 
the aircraft, and let (|)y(T,tlaf) and $y(f,tlaf) denote the con- 
ditional Instantaneous autocorrelation function and conditional 
Instantaneous spectrum of the response. Then, in Mark [ P ] it is 
shown that the conditional instantaneous autocorrelation function 
and spectrum of the response process are related to the correspond- 
ing characterizations of the input process w(t) and aircraft by 


<J)y(T,tla^) = 


J 

— 00 — 00 




<)>^(T-5,t-ula^)d5du 


(3.27) 


and 


$y(f ,t lOf.) 


00 

0^(f,u) $^(f,t-u|a^)du 

■ CD 


( 3 . 28 ) 


Substitution of Eq. (3.8) Into Eq . (3.28) yields an expression 
for the conditional Instantaneous response spectrum in terms of 
the spectra of the "slow" and "fast" components — i.e., Wg(t) and 
Wf(t) — of our turbulence model.: 


$y(f,t| 0 ^) = 


$^(f,u) (f)du + 


w 


$j^(f,u) 4>^^(f ,t-u]a^)du 


<!> (f) |H(f) I ' + 

w ^ 


$^(f,u) $^^(f ,t-ula^)du , 


(3.29) 


26 



where, in going to the second line, we have used Eq. lA on p . 26 
of Ref. 7 — applied to the deterministic function h(t) — l.e.. 


»oo 

$^(f,t)dt = lH(f) I"' , (3.30) 

_co 


whe re 


H(f) = 


h(t)e-i2^^^ dt 


(3.31) 


is the aircraft frequency-response function. 

Equation (3.29) is the desired spectrum input-response 
relationship. Notice that the first term in the right-hand side 
is the contribution to the response spectrum from the "slow" 
turbulence component Wg(t) and that this contribution has the 
usual form of a response spectrum. 


Series Expansion of Conditional Instantaneous 
Spectrum of Aircraft Response 


By substituting Eq . (3.12) into Eq . (3.29), and then per- 
forming the indicated integration term by term, we obtain a 
series expansion for the conditional instantaneous response 
spectrum $y(f,tlof): 




(f) |H(f)12 


N $^‘^^(f) 

+ y i — 

n=0 


O, (f,u) a (t-u)du 
h ’ n 


$^(f,u) (f,t-u)du , 


(3.32) 


where the third line in the right-hand side of the above ex- 
pression is the contribution of the remainder term. 
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The requirement of Eq. (3.20) Is that the fluctuations of 
a|*(t) occur slowly In comparison with the fluctuations of the 
von Karman component z(t). This locally stationary requirement 
Involves properties only of the turbulence. Obtaining our 
locally stationary representation of the aircraft response 
spectrum $y(f,t|af) further requires that fluctuations In a?.(t) 
be typically negligible over time Intervals comparable to the 
nominal duration of the aircraft Impulse response function h(t). 

We shall now derive explicit conditions that must be satisfied 
for a locally stationary response approximation to be valid. 

First, we recall from Eqs. (3.15) to (3.18) that the terms 
an(t) In Eq. (3.32) Involve only the time behavior of the tlmeT 
varying standard deviation Uf(t) of the ’’fast" turbulence com- 
ponent Wf(t) — see, e.g., Eq. (2.^). If Of(t) varies slowly In 
comparison with the nominal duration of h(t), then a few terms 
In the Taylor’s series expansion of an(t-u) about the Instant 
u = 0 should provide a good approximation to aj^(t-u) In Eq. (3-32) 

over the Interval of u In the Integral $j^(f,u) an(t-u)du, 

where |$h(f,u)| Is not negligible. This Taylor’s series expansion 
of an(t-u) Is 


a^(t-u) = a^(t) - u a^(t) - ••• 


= ? (t)u^ + if^^-,(t,u) , 


k=0 


‘M+1 


where 


a^^^t) 

n 




d^a^(t) 

dt^ 

a^(t) 


(3.33a) 

(3.33b) 

(3.3^a) 

(3.3^b) 


and where the magnitude of the remainder term satisfies 

(M+1) l,,|M+l 

I^M+1 I C (M+1) I ' (3.35) 

Here, 1 ^ ( C ) 1 designates the maximum value of the magnitude 

of a^^^^^(C) when ? Is varied over the range of values (t-u)£C£ t 
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or t£^£(t-u), depending on whether u is positive or negative, 
respectively. We may now substitute Eq . (3.33b) into Eq . (3.32) 
and Integrate, term by term. Defining 


m 


(k) 

h 



>00 

u^$j.^(f,u)du 

_oo 


(3.36) 


we find, upon carrying out this substitution and the resulting 
integration. 


<!>„ (f) |H(f) 

N 

+ 1 
n=0 


n! 

0 

II 

c 

+ 

$(n)(f) 

nl 

fOO 



2 




k=0 


$^(f,u) (t ,u)du 


$j^(f,u) Rj^^^(f,t-u)du , 


(3.37) 


where the last two lines are remainder terms that are of academic 
Interest only insofar as the present study is concerned. Equation 
( 3 . 37 ) is the desired series representation of the conditional 
instantaneous aircraft response spectrum. 


Requirements for Validity of Locally Stationary 
Aircraft Response Approximation 


The purpose of the present study is to determine conditions 
under which all terms in Eq. (3- 37) are negligible except for 
those corresponding to the locally stationary response approxi- 
mation. To accomplish this goal, we write out the double summa- 
tion in the second line in the right-hand side of Eq. (3-37) and 
make use of Eqs. (3.16) and (3.17). These operations yield 
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I ■ I 


^ = $ (f) |H(f)r 

^ s 

+ $2^^) lH(f) |2 - ^ a^Ct) + I 

dt 

+ I a^Ct) |H(f)|" - ^ a2(t)m^^^f) + | -^ a^ (t )m^^ ^ f )— • J 

d, t 

+ ••• , (3.38) 


where the first line on the right-hand side Is the contribution to 
the response spectrum $y(f,t|af) from the "slow" turbulence com- 
ponent Ws(t), the second line Is the expansion over k = 0,lj2j... 
of the term n = 0 In the double summation In Eq. (3-37) > the third 
line Is the expansion over k = 0,1,2,... of the term n = 2 In the 
double summation In Eq. (3.37)» and the +... In the last line 
represents terms n > 2 In the summation over n together with the 
remainder terms In the last two lines of Eq. (3-37) ■ Also, In 
writing out Eq. (3.38), we have used the relationship 

mj^°^f) = lH(f)l^ , (3.39) 


which is a direct consequence of Eqs . (3.30) and (3.36). We shall 
discuss Interpretation of the quantities m^^^(f) later In this 
report . 

The first line In the right-hand side of Eq. (3.38) requires 
no further discussion. The second line, which contains terms 
corresponding to n = 0 from Eq. (3.37), is the contribution to 
$y(f,t|of) from the locally stationary approximation to the con- 
ditional turbulence spectrum $wf (^>6 | erf ) . This Interpretation 

may be seen from Eqs. (3.12) and (3.19), where we remind the 
reader that the Index n In Eqs. (3.12) and (3.37) plays the same 
role. The third line In the right-hand side of Eq . (3-38), which 
contains terms corresponding to n = 2 from Eq. (3.37), contains 
the contributions to $y(f,t|af) from the first correction term 
to the locally stationary approximation to the conditional tur- 
bulence spectrum $i^^(f , 1 1 Cf ) . Furthermore, the first term within 

each of the brackets In the second and third lines In Eq . (3.38) 
provides the locally stationary approximation to the conditional 
response spectrum, whereas the second and third terms within each 
of the brackets In the second and third lines In Eq. (3.38) pro- 
vide correction terms to the locally stationary approximation to 
the conditional response spectrum. These second and third terms 
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within the brackets result from thie terms 
expansion of a^Ct-u) given by Eq . (3.33). 

to 


k = 1 and 2 in the 
Thus, the locally 

stationary response approximation to $y(f,t|of) that results from 
the locally stationary approximation to the turbulence spectrum 
lof) Is 




$y(f,t|oj.) « [ 4.^ (f) + a|.(t) 

s 


(3.40) 


which, of course, could have been written directly from the tur- 
bulence model of Eq. (3.1). 

Equation (3.38), together with the above comments, provides 
us with the criteria that must be met for the approximation of 
Eq. (3.40) to be valid. For the locally stationary approximation 
to the turbulence spectrum 4>^^(f,t|af) to be valid, the first 
term in the third line of Eq. (3.38) must be small in comparison 
with the first term in the second line. When Eq. (3.18b) is 
introduced into Eq. (3.38), this locally stationary excitation 
condition becomes 


d^ ilna^(t ) 
dt^ 


<< 32it^ 


$^(f) |H(f)|^ 
4>^^^(f)| |H(f)|^ 


(3.41) 


which must be satisfied in regions where af(t) is not negligible. 
Equation (3.41) is identical with Eq . (3.20), as expected, except 
for the terms lH(f)|^ in the right-hand side of Eq . (3.41). 
Inclusion of the terms lH(f)|^ relaxes the requirement of Eq . 
(3.20) to the extent that the right-hand side of Eq . (3.20) need 
be large in comparison with the left-hand side only over the 
range of frequencies where |H(f)|^ is not negligible. 

When the approximation of Eq . (3.41) is satisfied, our 

locally stationary response condition is satisfied if the second 
and third terms within the brackets in the second line of Eq . 
(3-38) are small in comparison with the first term. Using the 
fact that 


d 

dt 


£na^(t ) 


dt 


'f' 


a;.(t ) 


a^(t ) 


(3.42) 


and using Eq. (3.39)> the first of our locally stationary response 
conditions becomes 
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<< 




$^(f) 


(3.43) 


Notice that the Conditions of Eqs . (3.^1) and (3.43) both Involve 
derivatives of £na^(t). 


To relate the second locally stationary response condition 
to the derivatives of 5,na^(t), we note that 


dt 


5.na^(t ) 



d 

dt 


a 


2 

f 


and 


£na^(t) 

dt^ ^ 






(3.44) 


hence, rearranging Eq. (3.44), we have 


dt^ - d 


a^(t) 


dt 


ii-na~(t) + 
2 I 


dt 


(3.45) 


Prom Eqs. (3.39) and (3.45), it Is evident that the condition 
that the second term within the brackets In the second line of 
Eq. ( 3 . 38 ) be negligible In comparison with the first term may 
be expressed as 


— — £na^(t) + 
dt" ^ 


dt 


^.na^(t ) 


$z(f) |mj^°^(f)| 

^z(f) |nih^^(f)| ’ 


(3.46) 


which Is the second of our locally stationary response conditions . 
When Eqs. (3.41)^ (3.43), and (3.46) are all satisfied, the 

locally stationary conditional instantaneous response spectrum 
approximation given by Eq. (3.40) is valid. 
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The physical significance of the conditions of Eqs. (3-^1) > 
(3.^3)> and (3.46) will be discussed In a later section of this 

report, where expressions for m^^^(f) and mh^^(f) will be given 
In terms of the magnitude and phase of the aircraft frequency 
response function H(f). Methods for evaluation of the required 
turbulence metrics to test for the satisfaction of the condi- 
tions of Eqs. (3.41), (3.43), and (3.46) also will be derived 
later in this report. 

The importance of the locally stationary turbulence and 
response approximations has been pointed out by Sidwell — e.g., 
pp. 19-22 and p. 4l of Ref. 5. 
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AIRCRAFT RESPONSE EXCEEDANCE RATES AND 
PROBABILITY DENSITY FUNCTIONS 


Gaussian Property of Response Process 
Conditioned on the o^(t) Process 


In this section, the locally stationary Instantaneous re- 
sponse spectrum approximation given by Eq. (3.^0) will be used to 
derive expressions for aircraft-response exceedance rates and 
probability density functions.' In evaluating the usefulness of 
these results, we recall that the unit impulse-response function 
h(t). Introduced in the last section, may be regarded as describ- 
ing the response of any particular part of the aircraft modeled 
as a linear constant parameter structure. For example, h(t) may 
be chosen to describe the stress response at a critical section 
in a wing spar. 

In the derivations to follow, we shall assume that the "slow" 
component Wg(t) and the component z(t) in the turbulence model of 
Eq. (2.3) are both stationary Gaussian processes with zero mean 
values. When Wg(t) and z(t) satisfy this stationary Gaussian 
as'^umption, the response process y(t)t aonditioned on the process 
Oj!’(t)j is a zero mean strictly Gaussian {generally nonstationary) 
process. To prove this, we first note that each sample function 
of y(t) can be expressed as 


where 


y(t)|aj. = yg(t) + y^(t)ja^ , 


Y3(t) 


.00 

h(x) w (t-T)dT 

•mCO 


9 


and 


yp (t) |a^ 


.CO 

h(x) w^(t-x) Icj. dx 

.00 


(^. 1 ) 


(4.2) 


(4.3) 


where the vertical bars followed by af denote the fact that, in 
the stochastic process interpretation of the above results, the 
sample functions ap(t) are to be regarded as known or specified 
but otherwise unrestricted functions of time. According to 
Eq. (4.1), {y(t)|af} is the sum of two (Independent) processes 
{yg(t)} and {yf(t)]af}, since Wg(t) and z(t) are assumed to be 
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Independent. By assumption, . {wg (t ) } is a Gaussian process. Since 
any linear transformation of a Gaussian process results in a 
Gaussian process, e.g., see Cramer [S], pp. 312,31?, it follows 
that {yg(t)} is a Gaussian process. Thus, since the sum of any 
number of independent Gaussian processes is itself Gaussian, 
e.g., Cramer [5], pp . 316-317, it follows that the process 
{y (t ) 1 Of} will be Gaussian (but generally nonstationary), if the 
process {yf(t)laf} is Gaussian. However, according to Eq. (^.3), 
■[yf(t)laf} is a linear transformation of {wp(t)|af}; thus, 
{y(t)|af} will be Gaussian If {wf(t)|af} is Gaussian. Further- 
more, according to Eq. (2.^1), when a^(t) is considered as a known 
or specified function, wp(t) is a (memoryless) linear transforma- 
tion of z(t); hence, the process {wf(t)lof} is Gaussian. There- 
fore, the conditional response process {ytt)|af} is a (generally 
nonstationary) Gaussian process. This result does not depend 
on any of the locally stationary requirements of Eqs . (3.41) ^ 

(3.43), or (3.46). 


In the following work, we shall require the first-order 
probability density functions of (y(t)} conditioned on the 
process Of(t). According to the above comments, this probability 
density p(y|af-) is Gaussian; i.e.,* 


p(y |cfp 



(ii.i*) 


where (aylap)'^ is the standard deviation of the response process 
{y(t)|afl given that Of(t) is specified; i.e., (ay|af-)^ is the 
square root of the conditional variance: 


= E{y^ 1 op 


(iJ.5) 


We also shall require the fourth moment of y(t), given that Of(t) 
is specified. Prom Eq . (4.^) and known properties of the Gaussian 

density function — e.g., pp . 220-221 of Parzen [20] — this fourth 
moment is 


E{yPop = 3(aJ|op . (4.6) 

*We have used both p(***|af) and p(***lo|-) to denote the fact 
that the value of the process Of is specified. These are 
equivalent concepts. This equivalence is used in other nota- 
tions as well. 
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General Expression for Aircraft Response Exceedance Rates 


Numerous references in the- aeronautical literature have 
dealt with the mean rate of exceedances of an aircraft-response 
variable y(t) past some specified threshold level. Here, we 
shall deal with the expected number of crossings witH positive 
slope N+(y) of the response past a specified level y. It was 
shown by Rice [15] on pp. 189-193 of the Wax edition that, for 
stationary processes ^ one has : , 


N^(y) = 


>00 

yp(y,y)dy 

0 


(4.7) 


where p(y,y) is, the Joint probability density function of y and 
its derivative y. A derivation of this result also is given on 
pp. 45-47 of Crandall and Mark [52], The result of Eq. (4,7) 
applies to nonGausslah as well as Gaussian processes. Using the 
conditional joint probability density function p(y,y|a^) of the 
aircraft response y and its derivative y, we may formally express 
Eq . ( 4 . 7 ) as 


N^(y) = 


p(y>y kf ) 


p(apda^ 


dy 


0 0 


•oo 

■ 


. . 

0 

0 

r 


Np 


y p(y,ylopdy 


N^(y|a|.) p(apda^ 


p(apda|. 


(4.8) 


where we have defined 


N^(y kp 


► 00 

y p(yjykpdy 

. 

0 


(4.9) 


which we may regard as the expected number of threshold crossings 
with positive slope of the aircraft response past the specified 
level y, given that apt) is specified. 
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Several comments about Eqs. (4.8) and (4.9) are In order 
at this juncture. First, p(y,^|o|«) designates the joint pro- 
bability density of the response y and its derivative given 
that the stochastic /unction a^(t) is specified. Therefore, ' 
P(y»^lcTf) is a function of time. Furthermore, because of this 
fact, the integration involved in the formal determination of 
P(y>^) from p(y,J|af>) - i.e. , 


p(y»y) 


.09 

p(y»y|op P(apda* 

0 


(4.10) 



— must actually be regarded as an infinite dimensional integration 
over an infinite dimensipnal space a?*(t i ) , af-(t i ) , a|*(t s ) , . . . 
as for all i = 1,2,3,... is made to shrink to zero. 

Howeveri when interpreted in this manner, the operations in Eqs. 
(4.8) and (4.9) are valid. In particular, it follows from the 
above derivation that Eq. (4.9) is valid even though the . condi- 
tional process {y(t)|of(t)} is not stationary. 


In order to evaluate p(y,ylcTf) for the purposes of this 
report, we shall want to assume that {y(t)|af(t)} is locally 
stationary. When {y(t)|af(t)} is stationary , with zero mean value, 
it is known that E{y ylof} = 0; l.e., y and y are uncorrelated. 

It is shown in Appendix A of this report that when the conditional 
instantaneous spectrum of the aircraft response satisfies Eq. 
(3.40) — i.e., when the conditions of Eqs. (3.4l), (3-43), and 
(3.46) are satisfied — the conditional correlation coefficient 


, , ECyylap 
Pyyl^f 

y(aj|ap(ai lap 


(4.11) 


is negligibly small. Hence, for cases where Eq. (3.40) is valid 
and where the turbulence processes {ws(t)> and {z(t)} in Eq. (2.3) 
are Gaussian with zero mean values, the joint conditional pro- 
bability density p(y,y|af-) is two-dimensional Gaussian in the 
uncorrelated variates y and y. Hence, when the locally stationary 
condition of Eq. (3.40) is satisfied, we may apply Rice's well 
known formula till, p. 193 of the Wax edition, to the evaluation 
of Eq. (4. 9) : 


N^(y |op 


1 

27r 


X0,t\o^) 

iy(0,t la^) 



e 



(4.12) 
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where <t>y(T,tlaf) is the conditional instantaneous autocorrelation 
function of the response, defined in a manner completely analogous 
to Eq. (3.2a), and where the differentiation indicated in the 
brackets is with respect to t. 


To evaluate the quantity within the brackets in Eq. (4.12), 
we first note from the Fourier mate to Eq. (3*5) that the condi- 
tional instantaneous autocorrelation function <J>y(T,tlaf) is the 
Inverse Fourier transform with respect to f of the conditional 
instantaneous spectrum $y(f,t|af); l.e.. 



• 00 
I 

— oo 


(4.13) 


Hence, setting t equal to zero in Eq. (4.13) yields 


4« (0,tla^) = 


4> (f,t laf)df 


= cfy(t) (a^ . (4.14) 

Furthermore, differentiation of Eq. (4.13) twice with respect to 
T yields, after setting t equal to zero in the resulting expression 



-(2tt)2 


• 00 

f=‘4>y(f,t|a^)df 

_CO 


= \o^ 


(4.15) 


where the second line is a well-known result — e.g., pp . 190-192 
of the Wax edition of Ref. 11, or pp . 47-48 of Ref. 12 — and 
where the validity of the second line depends on the locally 
stationary property of Eq. (3-40). 

We may now use Eq. (3.40) to evaluate the right-hand sides 
of Eqs. (4.l4) and (4.15). Combining Eqs. (3.40) and (4.14) 
yields 
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(^. 16 ) 


aj(t) |aj = + a"(t) 


where 


j 

— 00 


(f) |H(f)lMf , 
s 


(i| 


and 


r2 § 


|H(f)lMf . 


(M 


Similarly, combining Eqs. (3.40) and (^1.15) yields 


(t) |a^= 0? + a^(t)al 

y ^ ^s ^ yz 

i 

where 

.oo 


a? = (2u)2 

f^ (f) 

w 

s 

.oo 

lH(f ) 1 Mf 

and 

,oo 


» I 

CM 

C\J 

II 

tsi 

CM 

D 

f^ ^^{r) 

.CO 

|H(f)|2df . 


(i» 




(4 


Finally, by combining Eqs, (4.14), (4.15), (4.l6), and (4.19) 
with Eq , (4.12), we obtain our final expression for the condi 
tional rate of exceedance crossings with positive dlope: 


N+(ykp = 


1 

27T 


(t )a 


yc 


aj(t)aj^ 




2[a§g+af (t)a^2^ 


(4 


which is valid whenever the locally stationary approximation 
Eq . (3-40) is valid. The quantities , and a? 

^s 


‘.17) 

. 18 ) 

.19) 

. 20 ) 

. 21 ) 

. 22 ) 

of 
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within the brackets in Eq. (4.22) are to be evaluated using 
Eqs. (4.17), (4.18), (4.20), and (4.21). When Eq. (4.22) is 
substituted into Eq. (4.8), we obtain a general expression for 
the mean exceedance rate of the aircraft response expressed in 
terms of the probability density function p(ai) of the stochastic 
function a|(t ) . 

Series Approximation of Aircraft Response Exceedance Rates 


Evaluation of ■ the expression for the response exceedance 
rate N+(y) given by Eq . (4.8) requires knowledge of the probability 
density function p(a|‘). However, for cases where the variance of 
Of — i . e . , 

(2) A 

^ E{[a|. - E(ap]2} 

= E{[apM - {E[ap}2 (4.23) 

— is small in comparison with the square of the mean E{a|>} of o^, 
a useful series approximation to N+(y) can be obtained, as we 
shall now show. 

Let us consider N+(y|af) as a function of the variable Of. 
Consider the Taylor's series expansion of N+(yIof-), where^the 
expansion is to be centered about the expected value of Of — l.e., 

= E{op . (4.24) 

We may formally write this expansion as 

N+(y|ap = (y lap (4.25a) 

= N^(ylap + Nj^^(y|ap (ap^) 

+ j Nj^^(y|ap (apap^ + ... , (4.25b) 

where we have defined 
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(i<.26) 




_ ^ dX(yl^p 


a^) 


d(ap 


k 


'°f = °f 


and 


N 


( 0 ) 


(y|ap = N^(yla^) 


(i*.27) 


Let us now substitute Eq. (^.25a) Into Eq . (^.8) and interchange 
the orders of Integration and summation: 


N^(y) 


00 


I 

k=0 


1 

k! 


N|*"^y|ap 


• 00 

(o^-op^ p(apd(op 

, 

0 


00 


= I 

k=0 


1 

k! 



(y kpy 



where 


y 





■°° 

(a^-ap^ P(apd(ap 

. 

0 


(4.28) 


(4.29) 


is the \cth central moment of the stochastic variable a?-. The 
reason for choosing the expansion center of the Taylor's series 

of Eq . (4.25) about the point may be seen from Eqs. 

(4.28) and (4.29). According to Eq.^(4.29), 

( 1 ) 

y„2 =0 ; (4.30) 

°f 


l.e., the term k = 1 In the expansion of Eq. (4.28) vanishes when 

2 

the center of the Taylor's series expansion is Of. Furthermore, 

It Is known — e.g.. Ref. 9, p. 175 — that the second moment of 

a|. Is minimized when taken about ap consequently, the term 

k=2 in the expansion of Eq. (4.28) Is minimized when the expansion 

center of the Taylor’s series expansion of Eq. (4.25) Is taken 

as ap The first few terms of Eq. (4.28) are 
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(^.31) 


N^(y) = N^(y|ap 


1 

2 


‘(2) 


N^^^Cylap + 


Discussion. It is immediately obvious using conditional 
expectations — e.g.. Ref. 8, pp. 55-56 — that when we consider 
of to be a stochastic variable, we may express the mean-square 
aircraft response E{y^} as 


E{yM 


r 


00 


E{y^ |ap p(apda^ 


0 


[a^ ] p(apda^ 

e ^ rr 


= a‘ 


y. 


2 2 


(^.32) 


where we have used Eq. (4.l6), with an obvious minor change in 
notation, and the definition of Eq . (4.2^). Consequently, It 
follows directly from Eqs . (4.22) and (4.32) that the first term 

series of Eq. (4. SI) is the response exceedance 
rate one would compute by assuming that the aircraft response 
y(t) is a stationary Gaussian process with a variance equal to 
the actual variance E{y^} given by Eq. (4.32). It follows that 
the terms k = 2,3,... in Eqs. (4.28) and (4.31) are correction 

terms to the Gaussian approximation E_^(y\a^) of the true air- 
craft exceedance rate N^(y)y where these correction terms take 
into account the nonGaussian nature of the aircraft response 
that is caused by the stochastic variations in the quantity 
Ojf(t) in the model of Eq. (2.3). 

The behavior of the first correction term. 




( 2 ) 


(y|a^) y 


( 2 ) 


to the Gaussian approximation N+(ylaf>) in Eq. (4.31) Is of 
special Interest, since this correction term governs the behavior 

( 2 ) 

of the deviation from Gaussian behavior N+(y) when y „2 is small 

— . ^ f 

In comparison with It is shown in Appendix B of this 
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can be expressed, exactly, as 


report that N^^^(y|a|.) 


N 


j^^(y|ap = 2 N^(y|ap (y |crp 


(2) 


(^. 33 ) 



and where we have used the notation 




= 

+ op^ 

(4.35a) 


y ' f 



and 

a|lo“ 


+ 0p| 

^ z 

(4.35b) 

When Eq. 

(4.33) 

is substituted into Eq. 

(4.31), and only one cor- 


rectlon term to the Gaussian approximation is retained, we obtain 
N^(y) s N^(y|ap [1 + Q^^^(ylap] ; (4.36) 

hence, we see that the quantity Q' '^(y|c7f) is of the nature of a 
multiplicative correction factor to the Gaussian approximation 

N+(y lop , 

Often — e.g.. Ref. 13 — the logarithm of the normalized 
exceedance rate N+(y)/N+(0) is plotted. Dividing Eq. (4.36) by 
N+(0) and taking the natural logarithm of the resulting expression 
gives 


43 



Jin 


N+(y) N+(y|ap 


N+(0) 


= An 


( 2 )^( 2 ) 


N, 


^-Qy- + An[l + y ^2 Q (y|ap] 


N^(y|ap (2) (2) — 

N^(o) * Woi « • 


( 2 ) ( 2 ) — 

for y ^2 Q (y|ap < 0.2 , 


(i<.37) 


where. In going to the second line, we have used the first term 
in the Maclaurin expansion of An(l+x) — i.e., An(l+x) » x, which 
has approximately ten percent error when x = 0.2. It follows from 


Eq. (4.Z7) thatj for 
function ^ (y\a^j.) 


sufficiently small values 
governs the deviation of a 


of y 
plot 


( 2 ) 



the 


inlN_^(y ) /N _^( 0 ) 2 from the quantity ln\_E ^(y\o^) /E _^( 0 )2 which is, 

except for a possible constant factor, the result predicted if 
one assumes that y(t) is stationary and Gaussian. 


In order to discuss further the result of Eq . (4.37), we 
shall restrict our attention to the case where the aircraft re- 
sponse to the component Wg(t) of our model of Eq. (2.3) is neg- 


ligible. This assumption implies that a 


y< 


0 In Eqs. (4.35a) 


and (4.35b). Consequently, for this llmitin 
from Eqs. (4.35a) and (4.35b) whe■^.a^ and 
zero, 


case, we obtain 
ai are set equal to 
hs 


a 




y. 


a? I 
y ' f 


(4. 38a,b) 


(4.38a,b) into Eq . (4.34) and writing Oy for 
Oylof, we obtain, after simplification. 


Substituting Eqs 

r2 I „2 




^ 8(ap" 


y-ly--h 

y 


(4.39) 


which is valid only when Oy = 0 and a\ = 0. Notice that, for 

2 ^ S V S 

this case, Oy = E{of}Oy ; hence, the quantity Oy in Eq. (4.39) is 
the true variance of tnl response variate y(t). Furthermore, 
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according to Eq. (^.39), we have Q^^^COlaf-) = 0. Hence, 
according to the approximation of Eq . (^.36), we have N+(0) = 

N+(0laf). It follows that, for the case under consideration, 
the first term in the right-hand side of Eq. (4.37) can be 


written as S.n[N+(y | of )/N+( 0 | ] which is, exactly, the loga- 

rithm of the normalized exceedance rate for a stationary Gaussian 
process. Using Eq. (4.22) to evaluate this term and writing 

Oy for ofOy^ (recall that Oy^ = 0), we have for the right-hand 


side 


of Eq. (4.37) , 


An 


N^(y) 

nTToT 



( 2 ) 

1 ^ 


ii 


y 



(4.40) 


where we have used Eq. (4.39) and where, according to Eq. (4.37), 
the correction term, which is the second term on the right-hand 
side of Eq. (4.40), should be accurate to within about ten percent 
whenever that term is less than about 0.2. 


Equation (4.40) is the desired result. If Jln[N^(y )/N_j.(0 ) ] 
is plotted on the ordinate versus y^/Oy on the abscissa, we see 
that the first term in the right-hand side of Eq. (4.40) is a 
straight line with slope of minus one-half and ordinate inter- 
section at (y^/Oy) = 0. This first term is the exact result for 
stationary Gaussian processes y(t). The second term on the 
right-hand side of Eq . (4.40) is a parabola in the quantity 
y^/Oy. This second term is zero at (y^/cy) = 0 and at (y^/dy) = 
Between these two values, this correction term is negative (since 


( 2 ) 

y 2 is necessarily positive or zero). 


For values of (y^/Oy) > 4 


4. 


9 


this correction term is positive. Consequently y for threshold 
levels between 0 < \y\ < SOy, Eq. (4.40) indicates that the first 

(stationary Gaussian) term overestimates the number of threshold 
crossings; whereas y for thresholds \y\ > 2Cy, Eq. (4.40) indicates 

that the first (stationary Gaussian) term underestimates the 
number of threshold crossings of processes y(t) with time-varying 
variance . * 


Equation (4.40) provides a theoretical prediction of the 
concave shape of aircraft-response logarithmic exceedance plots 
for the turbulence model described by Eq. (2.3) for the case 
where Wg(t) = 0. The "strength" of the parabolic correction 
(in the variable y^/Cy) to the stationary Gaussian response 
result, given by the first term on the right-hand side, is 
governed by the coefficient 


*For y = 0, 
stationary 

Q^^^(Ol^) 


it is evident from Eqs . (4.36) and (4.39) that the 
Gaussian term gives the correct value for N+(0) since 

is zero. 
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(ii.illa) 


u^2) 

E{[a^-E(0p]^} 

{ECoJ.]}" 

_ E{[a|]n - {E[a|]}2 

E{[ap}" 


(4.lllb) 


which Is the square of the ooeffioient of variation (e.g.. Ref. 9) 
of the time-varying variance af(t) in the turbulence model of 
Eq. (2.3). 


Aircraft Response Probability Density Functions 


It is evident from Eq. (^.12) that, when the aircraft 
response y(t) is a stationary Gaussian process, the mean rate 
of "up crossings" N+(y) of the level y is proportional to the 
probability density function (pdf) p(y) of y(t). We shall show 
in this section that the proportionality between N+(y) and p(y) 

Is no longer maintained when the turbulence excitation Is modeled 
by Eq. (2.3). 


First, we shall derive a series expansion for the first 
order pdf p(y) of the response process {y(t)} using a method 
completely analogous to that used to obtain our series representa- 
tion of N+(y), Denoting by p(ylof) the conditional pdf of y 
given that is specified, we have 


P(y) 


.00 

p(ykp p(0|)da| 

0 


(^.i»2) 


We now formally expand p(y|0p) in a Taylor's series in the 
variable 0 ^ about the expansion center E{ 0 ^} = 0 ^: 


p(ylop = I — p^^^(y|0i) (0f-0p^ , 

^ U-=0 k! Ill 


(^.^3) 


k=0 k! 


where we have defined 


iiS 



(H.iiin 


and 




P<“’(y|pp = p(y|o5.)| . — . 


— ^ 4 <aS(yl°p 

cJ(pp’' 


_ _2 
^ f ^ f* 


2 2 


(^.^5) 


Next, we substitute Eq. (^.43) Into Eq. (4.42) and Interchange 
the orders of integration and summation, finding thereby that 


p(y ) 


= y — 

l k ' 
k =0 


p“'’(y 


op 


(k) 




(4.46) 


where we have used the definition of Eq. (4.29). The motivation 

for using Of as the expansion center of the Taylor’s series ex- 
pansion of Eq. (4.43) is the same as that described earlier in 
the text between Eqs. (4.30) and (4.31). We may write out the 
first few terms of Eq. (4.46) as 

p(y) = p(y|crp + j p^^^(y|cTp + ••• , (4. 47 ) 

since the term in Eq. (4.46) corresponding to k = 1 is Identically 
zero because we have chosen our Taylor's series expansion center 

at = a^. 

The first term in the right-hand side of Eq. (4.47) is the 
Gaussian approximation to the pdf of y(t) given by Eq . (4.4) with 

^f “ according to Eq . (4.35a), we have 


l^f = 


+ ola^ 


f y. 


E{y^} = a 


(4.48) 


The second term in the right-hand side of Eq. (4.47) is the 
"first-order" correction term to the Gaussian approximation 
provided by the first term. When the coefficient of variation 

lv^ 2 ^/(op^i^ of the random variable ai is sufficiently small, 
0^1 I 

this first-order correction term will exhibit the principal devia- 
tion of p(y) from the Gaussian approximation p(yja^). 
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( o\ 

It Is shown in Appendix C of this report that p^ 
can be expressed, exactly, as ^ 


P^^^(y|ap = 2 p(y|a|) U^^^(y|ap , 


(^.49) 


where 


(aj 

u^^^ylap - I 


1 

F 


(aj - 

z 


1 - 


2y- 


+ ^ 




- 1 




"yl^f 


- 6 


°yl°f 


+ 3 


(^.50) 


Substitution of Eq. (^.49) Into Eq. (4.47) yields the two-term 
approximation to p(y): 


p(y) 


p(y 


af) 


1 + 


Of 


0 f) 


(4.51) 


Disoussion. It Is in the nature of the derivation of Eqs . 
(4.46), (4.4>7), and (4.51) that, when y^f^ Is sufficiently small, 

O £• 

the right-hand side of Eq. (4.51) should provide a good approxi- 
mation to p(y). It Is of considerable Interest that, using an 
apparently completely different line of reasoning — namely, the 
Gram-Charller series — e.g., Cramer, Ref. 9 — one obtains pre- 
cisely the same correction to the Gaussian approximation given 
by Eq . (4.51). This equivalence Is most readily seen by a 
comparison of Eqs. (4.50) and (4.51) above with Eq. (8-111) on 
p. 272 of Papoulls [74], where we note that the quantity within 
the brackets in Eq . (4.50) above is the Hermlte polynomial of 

degree four In the variable y/(ayla^)'^. 


To illustrate the equivalence of Eq. (4.51) to the Gram- 
Charller series with the same number of terms, we must show, by 
comparing Eqs. (4.50) and (4.51) above with Eq . (8-111) of 
Ref. 14, or with Eqs. (17.6.5) and (17.6.6) of Ref. 9, that 

the coefficient of* excess Yy of y — i.e.. 


Y 


( 2 ) 

y 



(4.52) 
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(where 



is the fourth central moment of y) — Is related to 


V 


( 2 ) 

Of 


by 


(2 






(^. 53 ) 


or, equivalently, that 


Y 


( 2 ) 

y 





( 4 . 54 ) 


The equality in Eq . (4.54) is proved in Appendix D. 

The main purpose of the above discussions of response ex- 
ceedance rates and probability density functions has been to 
Illustrate jthe importance of the pdf p(af-) of the -^ime-varylng 
variance Of(t). The fundamental Importance of p(Of) in aircraft- 
response statistics is clearly illustrated by Eqs, (4.8) and 
(4.42). Furthermore, it is evident from the series expansions 

of Eqs. (4.28) and (4.46) that the (central) moments 

Of 

k = 2,3,4,... constitute th| most important set of parametria 
descriptors of the pdf of Cf. Methods for computing these 
quantities from turbulence records will be described in. Sec. 6 
of this report . 

Use of conditional probabilities in computing an expression 
for exceedance rates analogous to Eq . (4.8) of this report, but 
using p(of) rather than p(af), has been used recently by Sidwell 
[IS]. 


Comparison of Exceedance Rates and Probability Density Functions 


It is evident from Eqs. (4.4), (4.12), (4.36), and (4.51) 
that, considered as a function of the response variable y, the 
Gaussian approximations to the exceedance rate N+(y) and the pdf 
p(y) provided by the first terms in the right-hand sides of 
Eqs. (4.36) and (4.51) are proportional to one another; l.e., 
each has the form C exp[-y^/(2a^ | a|.) ] , where C is a constant. 
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However, this same statement* Is not true of the approximations 
to N+(y) and p(y) provided by the right-hand sides of Eqs. 
( 4 . 36 ) and (4.51), respectively. In particular, for the case 
where = 0, it follows from Eqs. (4.36) and (4.39) that the 

zeros of the correction term “^(ylaf) occur at (y^/Cy) = 0,4; 
whereas, from Eqs. (4.50) and (4.51) it follows that the zeros 

of the correction term U^^^(ylaf) occur at (y^/oy) = (3±/U) » 
0 . 55 , 5 . 45 . Thus, differences between the shapes of observed 
exceedance functions N+(y) and probability densities p(y) are 
an indication of nonGausslan behavior of the turbulence excita- 
tion process w(t).* 


*0n p. 34 of Ref. 5, Sidwell has pointed out that N+(y) and p(y) 
will, in general, have different functional forms. 
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LOCALLY STATIONARY RESPONSE REQUIREMENT INTERPRETATIONS 
AND RELATIONSHIPS TO TURBULENCE MEASUREMENTS 


Interpretation of Locally Stationary Instantaneous 
Response-Spectrum Requirements 


Three requirements for the validity of the approximation of 
Eq . (3-^0) for computation of the conditional instantaneous 
aircraft response spectrum were arrived at in Sec. 3 of this 
report. These requirements are given by Eqs. (3-^1), (3.^3), 
and (3.^6). 

The first requirement, given by Eq. (3-41), has already been 
discussed in Ref. 6 in terms of a spatial variable x = Vt , where 
V is the aircraft speed. [See Eq. (^.82) of Ref. 5.] When the 
spectrum of z(t) of the model of Eq . (2.3) is the von Karman 
transverse spectrum with integral scale Lz, the condition of 
Eq. (3.^1) can be written as 

d^ilnai(t) ..2 

< 0.08 — , (5.1) 

dt^ Lj 

where we note that ina^ = 2Zna. The condition of Eq. (5.1) means 
that the second derivative of Znof(t), measured on the time scale 
t' E L 2 /V, must be less than or equal to O.O 8 . This requirement 
means, roughly, that the relative changes in the time-varying 
variance a\(t) of the "fast turbulence component" Wj^(t) in Eq. 

(2.3) must'' occur slowly when measured on the time scale t' = L^/V 
defined by the aircraft speed V and the integral scale Lz (i.e., 
nominal correlation interval) of the component z. We refer the 
reader to the discussion of Eq . (H.82) in Ref. 6 for further 

interpretation of Eq. (5.1). 

Equation (5.1) is independent of properties of the aircraft 
unit-impulse response function h(t), while the second and third 
requirements are dependent on properties of h(t). The second re- 
quirement given by Eq . ( 3 . ^3 ) , depends on the magnitude of the ratio 

(f )/m^^ ^ (f ) , where these two quantities are defined by Eq. 

(3-36). The quantities m^^^(f), k = 0,1,2,... were defined 

earlier in Ref. 16 by the author of this report for use in 
another context. It is shown on p. 282 of Ref. I 6 that the 

quantity mj^^^ (f )/m^^^ (f ) is the group delay T^(f) of the aircraft 

unit-impulse response h(t); i.e.. 
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(5.2) 


T^(f) 


A 


(1 


To 


(f) 

(f) 


If we express the Fourier transform, Eq. (3.31) of h(t) in terms 
of Its magnitude and phase 


i0^(f) 

H(f) = |H(f)| e ^ 


(5.3) 


then T^(f) can be expressed in terms of the phase by 


Th(f) = 




2ir 


df 


( 5 .^) 


The group delay Tj^(f) of h(t) can be interpreted physically as 
an energy-spectrum-weighted frequency-decomposition of the time 

centroid t^ of h^(t); l.e., 

th^ (t )dt 


h^ (t )dt 


See Eq. (165) of Ref. l6 and the accompanying discussion. 

Substituting Eq. (5.2) into Eq. (3.^3)j the second of our 
locally stationary response spectrum conditions becomes 

<< , (5.6) 


dt 


Zn ai(t) 



where we have Ignored the ratio 4>2 (^')/^z ^ (3.^3). The 

requirement of Eq. (5.6) means, roughlj^, that the relative 
changes in the time-varying variance a^(t) of the "fast turbu- 
lence component" Wjf(t) in Eq. (2.3) must be small in comparison 
with unity when measured over time intervals of the order of 
the group delay '^^(f) of h(t). 


For unit-impulse response functions h(t) that are even 
functions of time - i.e., h(-t) = h(t) - it can be shown from 
Eq. (5.4) that Th(f) = 0. [Recall the cho1 ce of our time origin 
of h(t) as defined by Eqs. (3.23) and (3.24).] For such 


52 



impulse-response functions, Eq. (5.6) is always satisfied; 
hence, the third condition given by Eq. (3.46) is required. 

This third condition depends on the ratio ^ ( f )/m]^^ ^ ( f ) , 

which is an energy-spectrum-weighted frequency-decomposition 
of the second moment of h®(t) normalized to unit area. Since, 
by our choice of the time origin of h(t), the time centroid of 

h^(t) is zero, it follows that m^^ ^ ( f )/m^^ ^ ( f ) is a frequency 

decomposition of the second central moment of h^(t) (normalized 
to unit area). See p. 286 of Ref. 16. It is shown on pp . 283 

and 284 of Ref. 16 that the quantity m^^ ^ (f )/m^*^ ^ (f ) can be 

expressed in terms of the magnitude and phase of H(f) by 


( 2 ) 


W) 


(f) 

(f) 


. 1 

d0^(f)' 

2 \ 

1 dnn|H(f ) 1 1 

4it^ ( 

df 

^ df^ j 


(5.7) 


Because m^^ ^ ( f ) /m^*^ ^ ( f ) is an energy-spectrum-weighted frequency- 

decomposition of the second central moment of h^(t) [with h^(t) 
normalized to unit area], the requirement of Eq. (2.46) means, 
rqughly that the relative changes in the time -varying variance 
Ojf(t) of the "fast turbulence component" wj^(t) in Eq. (2.2) must 
occur slowly when measured over time intervals comparable to the 
nominal duration of h^(t). Since the nominal duration of h^(t) 
governs the "correlation interval" of the response process, when 
the requirements of Eqs . (3.41), ^3.43)j and (3.46) are satis- 
fied, the fractional changes in Of(t) should be small over in- 
tervals comparable to the "correlation interval" of the response 
process {y(t)}. 


Locally Stationary Requirements Expressed in Terms 
of Mean-Square Response 


The left-hand side of each of our locally stationary re- 
quirements of Eqs. (3.41), ( 3 . 43 ), and (3.46) depends only on the 
behavior of £nof(t). Furthermore, when like terms are cancelled 
in the right-hand sides of these same equations, it is evident 
that the right-hand side of Eq. (3.41) depends only on the 
spectrum of the turbulence component z(t), whereas the right- 
hand sides of Eqs. (3-43) and (3-46) depend only on properties 

m^^^(f), and m^^^(f) of the aircraft impulse-response 
function h(t ) . 
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By Integrating the serTes expansion of Eq. (3-38) over 
_oo < f < ooj we obtain an expression for the conditional mean- 
square aircraft response E{y^(t)|af}. It Is Immediately evi- 
dent that, after this Integration Is carried out, the three 
locally stationary requirements of Eqs . (3.^1), (3*43), and 
(3.^16) are replaced by 


and 


d^ ilna^Ct ) 


df 


I t^(f) |H(f) l^df 


<< 32 ir‘ 


• 0 

. 


(f ) |H(f) I Mf 


(5.8) 


dt 


»oo 

$^(f) mj^°'’(f)df 


<< 


• 0 


^^(f) m^^^(f)df 


(5.9) 


dt 


- £na^(t) + 
2 f 


^ £na^(t) 


<< 2 


j 

.00 


‘I'zXf’) m^°^(f)df 


$z(f) mj^^^(f)df 


(5.10) 


The right-hand side of each of the above three conditions now 
depends both on the spectrum $ 2 (^) turbulence component 

z(t) and on one or more of the quantities m^^^(f) = |H(f)|^, 

f 1 ) (2) “ 

m^ '^(f), and m^ "^(f) that are determined from the (complex) 

aircraft frequency-response function H(f). However, the right- 
hand sides of Eqs. (5.8), (5.9), and (5.10) have the advantage 
that, for a given turbulence component z(t) and aircraft, each 
Is defined by a single number only. These conditions are some- 
what less restrictive than those defined by Eqs. (3.^1), (3.^3)j 
and (3.^6). 
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Expression of Requirements in Terms of Autocorrelation 

Function of Jlna^(t) 

The left-hand side of the three requirements of Eqs. (3.^1), 
(3.^3)s and (3.^6) or Eqs. (^.8), (5.9), and (5.10) are expressed 
In terms of the function £naf(t), which, until now, has been 
considered to be known or specified. However, In the turbulence 
model of Eq. (2.3), cJf(t) Is assumed to be a sample function 
drawn from a stationary stochastic process. We shall now ex- 
press the left-hand side of the above three locally stationary 
conditions In terms of the autocorrelation function of )lnaf(t). 

For brevity of notation, let us define 

v(t) = Jlna^(t) . (5.11) 


Denote the autocorrelation function of v(t) = ilna|*(t) by Rv(x); 
1 . e . , 


R^(t) = E{v(t) v(t+T)} . (5.12) 

Then, denoting by a prime the derivative of a function with 
respect to Its argument, we have 

H^(t) = E(v(t) V'(t+T)} 

= E{v(t-T) v'(t)} , (5.13) 

since (v(t)} Is a stationary process. A second differentiation 
yields 


R^(t) = E{-v'(t-x) v'(t)} 
= -E{v' (t) V' (t+T)} 


(5.1^) 


Hence , 


-r;;( 0) = E{[v’(t)]M ; ( 5 . 15 ) 

l.e., the mean square value of v'(t) Is equal to -Ry(0). Further- 
more, replacing v(t) by v’(t) In Eq. (5.12) yields, using the 


I I III I ^ra I ■■ 


I II I III 


result of Eq. (5.15), 

= E{[v»(t)]M , (5.16) 

where the superscript on the left-hand side represents the fourth 
derivative. The results of Eqs. (5.15) and (5.16) are, of course, 
well known — e.g., p. 21 of Ref. 17. Since v(t) = ^.na|*(2) is a 
stationary process, v'(t) and v"(t) are stationary processes with 
zero mean values. Hence, the standard deviation of v'(t) is 

[-Rv(O)]'^, and the standard deviation of v”(t) is [Ry^^(O)]^. 

Suppose now, by our << signs in Eqs. (5.8), (5.9), and (5.10), 
we mean that the first neglected terms in Eqs. (3-38) should 
be a factor of about three less than the terms they are being 
compared with when the stoahastia left-hand sides of the three 
conditions of Eqs. (5.8), (5.9), and (5.10) are statistically 
rather large, say at a "level" of three times their standard 
deviations. These two factors of three give us a composite 
factor of three times three, which we shall round off to ten. 
Consequently, the signs << in Eqs. (5.8) to (5.10) require a 
factor of about ten for their validity, when the left-hand sides 
are replaced by their standard deviations. For engineering 
applications, choice of a factor of ten as this minimum value 
for << is about right — l.e., not overly conservative. Hence, 
using Eqs. (5.16) and (5.15), we may express the requirements 
of Eqs. (5.8) and (5.9) as 


•oo 

$^(f) lH(f) I ^df 

[R^^^O)]^^ < 3.27t2 

(f ) |H(f ) I ^df 
z 


and 


[-r;;(o)]''^ 1 ^ 


w / \ 

$^(f) m^°^(f)df 

OO 

j mj^^^f)df 

«co 


(5.17) 


(5.18) 


With the condition of Eq. (5.10), the situation is some- 
what more complicated. The expected value of the left-hand side 
(without absolute value signs) is not zero in this case. 
Furthermore, we require the standard deviation of the sum of the 
two quantities on the left-hand side, and these two quantities 
are not statistically Independent. 
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To provide an expression for the standard deviation of 
the left-hand side of Eq. (5.10) that is of practical use, it 
is necessary to assume that v'(t) = d(Ana^)/dt is a Gaussian 
process. This assumption would not normally be Justified; 
however, for the kinds of Inequalities required here, the 
assumption should not cause appreciable error. When this 
Gaussian assumption for v'(t) is made, it is shown in Appendix 
E of this report that the required variance for the left-hand 
side of Eq. (5.10) is 

Var{v"(t) + [v’(t)]='} = R^^^(O) + 2[R»(0)]^ . (5.19) 

Furthermore, the mean value of the left-hand side of Eq. (5.10) 
(with absolute value signs removed) is 

E{v"(t) + [v'(t)]M = E{[v'(t)]M = -R»(0) , (5.20) 


according to Eqs . (5.11) and (5.15). Hence, using the three 
times the standard deviation criterion discussed earlier, and 
recognizing that this must be added to the mean value, we shall 
want for the left-hand side of Eq . (5.10) 

{2[R;;(0)]2 + R^^^O)}"'® - I R^(0) 



which is 1/3 times (mean plus three times standard deviation) 
and where we have recognized the fact that R(J-(0) is negative. 
Using the "three times three equals ten" rule discussed earlier, 
the requirement of Eq . (5.10) can now be expressed as 


1 

3 


r;;(o) 




+ 


R^^^O) j 

[r;;(o)]2 j 



1 

5 


«>z(f) m^°^f)df 
— 00 

^z^^^ ni]^^^(f)df 

— OO 


9 


(5.22) 

which is our third requirement. Equations (5.17), (5.18), and 

(5.22) are the conditions required for confident engineering use 
of the locally stationary approximation of Eq. (3.40) to compute 
the mean square aircraft response as was done, for example, in 
Eqs. (4.16), (4.17), and (4.18). 
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Method for Evaluation of Autocorrelation 
Function of £na^(t) 

The left-hand sides of Eqs. (5.17), (5.18), and (5.22) 
depend only on the second and fourth derivatives of the auto- 
correlation function Rv(t) of v(t) = £naf(t) evaluated at t = 0. 
Hence, to determine If these three conditions are satisfied, 
we must evaluate the autocorrelation function Rv(t) from measured 
turbulence records. This problem will be treated now. 

Let us return to the turbulence model of Eq. (2.3). Con- 
sider a high-pass filtered version w^(t) of w(t), where the 
filter cutoff frequency is sufficiently high to remove, for 
practical purposes, the "slow process" Ws(t). Hence, we may 
write this high-pass filtered version of w(t) as 

w^(t) = Aa^(t) z^(t) , (5.23) 

where Zh(t) is a high-pass filtered version of the component z(t) 
in Eq. (2.3) assumed normalized so that E{zh) = 1, A is the con- 
stant required for this normalization, and where it has been 
assumed in writing Eq . (5.23) that fluctuations in Cf(t) occur 
slowly in comparison with fluctuations in z(t). It is presumed 
that w^(t) would be generated from a turbulence recording by 
digital filtering. Thus, we shall assume that the record w^(t) 
is available. 

Let us now square Eq. (5.23) and then take its natural 
logarithm. These operations yield 

£n 'w^(t) = iln A^ + £n a^(t) + £n z^(t) . (5.2^) 

By assumption in our turbulence model of Eq. (2.3), (cf(t)} and 
{z(t)} are independent processes; therefore, the same is true of 
{£n a|'(t)} and {£n z^(t)}. Hence, we may immediately write the 
covariance of Eq. (5.2^), taken at observation times t and t + t, 
as 

cov[£n w^(t), £n w^(t+r)] = cov[£n a^(t), £n a^(t+T)] 

+ cov[£n z^(t), £n z^(t+x)] , 

(5.25) 
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where we have used the fact that the covariance of the sum of 
two independent variables is equal to the sum of the covariances 
and where we note that the dependence on the constant A has 
disappeared. According to the material contained in the pre- 
vious section, our interest is in the autocorrelation function 
v(t) = Jin af(t), which is related to the covariance function of 

v(t) = «.n o|.(t) by 

R^(t) = cov[v(t), v(t+T)] + {E[v]}^ . (5.26) 

Hence, we wish to solve Eq. (5.25) for cov[Jln a|.(t),Jln a|>(t+T)]; 

i . e . , 


cov[iln a^(t), Jin a^(t+x)] = cov[Jln w^(t), Jin w^(t+x)] 

- cov[Jln z^(t), Jin z^(t+x)] 

(5.27) 

The first term on the right-hand side of Eq. (5.2?) can be 
numerically computed directly from the high-pass filtered version 
Wj^(t) of the original sample w(t). Thus, in order to compute 
Rv("r)j V = Jin Of, we need an expression for the second term on 
the right-hand side of Eq . (5.27). To obtain this expression, 
we first recall that, by assumption, z(t) is a stationary 
Gaussian process; therefore, z^(t) also is a stationary Gaussian 
process. In Sec. 6 of this report, we shall describe a method 
for computing the power spectral density of z(t) from the mea- 
sured turbulence samples w(t). Let us denote this power spectral 
density by <J> 2 (f), and denote by Hh(f) the high-pass filter 
complex frequency-response function. Then, when the gain of 
this filter is appropriately adjusted* so that E{z^} = 1, the 
power spectral density Ozv,(f’) of ZK(t) can be computed from 
^^(f) by 

<I>2 (f) = 'J>2(f) r • (5.28) 


Hence, from a measured turbulence record w(t), we shall have at 
our disposal the means to compute the power spectral density 
^z^(f) of the component zh(t). Therefore, by using Eqs. (5.26), 


*It will become evident later that this assumption need not be 
satisfied in practice. 
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(5.27), and (5.28), we shall have at our disposal the means to 
compute Rv(t) If we can obtain cov[iln zA(t), An z^(t+t)] from 
the power spectral density $zj^(f) of Zh(t). 

Let us define 

u(t) ^ iln z*(t) . (5.29) 


Then, the covariance of u(t), u(t+r) Is related to the autocor- 
relation function R^^(T) by 

covCu(t), u(t+x)] = R^(t) - {E[u]}^ . (5-30) 

Thus, If we can compute the autocorrelation function Ru(t) of 
the function u(t) = *.n z^(t), we have at our disposal the means 
to compute R^(t). 

The problem, therefore, is to compute the autocorrelation 
function Ru(t) of the natural logarithm of z^(t), given that we 
know the power spectrum ^>zj^(f) of {zh(t))and under the assumption 

that {zh(t)} is a stationary Gaussian process with zero mean and 
unit variance. This problem Is solved in Appendix P, where It 
Is determined that 

R (t) = 2 arcsin^ p (t) + (C + Jin 2)^ , (5. 3D 

U 


where C Is Euler's constant — l.e., 

c = 0.577215665... 


(5.32) 


— and where Pzj^^'*^) correlation coefficient* of the process 

z, (t); l.e.. 


P., (t) = 


R^ (T) 
^h 

Rz (0) 
^h 


(5.33) 


*The formulation of Eq. (5.33) negates the requirement that the 
high-pass filter gain be adjusted so that E{z^} = 1. 
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From Eqs. (5.31) and (5. 33) 3 .. it immedlatelj’ follows that 

cov[£n zMt), In z^(t+x)] = 2 arcsln^tR (x)/R ( 0 )]. ( 5 . 3 ^) 

h h 


Finally, we note that Rzj^('’^) computed from ^>z^(f) by 


R„ (x) 
^h 


$ (f) df 

^h 


<I>z(f) l^h^f) I"" df 


(5.35) 


according to Eq. (5.28). Equations ( 5 . 26)3 ( 5 . 27)3 ( 5 . 34)3 and 

(5.35) ooltectively provide the means to compute the autocor- 
relation function Fy ( x>) of v(t) = in o^(t) from the high-pass 
filtered version W}i(t) of the turbulence record w(t) and the 
power spectral density ^^(f) of the turbulence component z(t): 

R^(x) = {E[iln 0 ^(t)]}^ + cov[J!.n w^(t), «-n w^(t+x)] 

- 2 arcsln^[R (x)/R (0)] , (5.36) 

where Rzj^(t) Is to be computed from ^z(f) and the high-pass 

filter frequency-response function by Eq . (5.35). Notice that, 

in computing the derivatives Ry(x) and R^^^(x), the constant term 
{E[£n Of(t)]}^ In Eq. (5-36) becomes irrelevant. Furthermore, 
Insofar as computation of the derivatives of Rv(t) is concerned, 
we may replace cov[Jln w^(t), in Wh(t+x)].by the autocorrelation 
function of 2.n w^(t) for this same reason. The computation of 
cov[Jln w^(t), Jin wn(t+x)] is to be carried out directly from 
the high-pass filtered version of the turbulence record w(t). 
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METHODS FOR COMPUTATION OF TURBULENCE METRICS FOR 
PREDICTION OF AIRCRAFT RESPONSE STATISTICS 


In addition to the turbulence metrics Rzj^(t) and Gov[iln w^(t), 

£n Wh(t+x)] used to compute the second and fourth^derivatives of 
the autocorrelation function Rv(t) of v(t) = Zn Of(t), which are 
required for verification of the locally stationary aircraft 
response assumptions, the material described in Sec. 4 of this 
report delineates several other turbulence metrics necessary to 
predict the response exceedance rates and probability density 
functions of an arbitrary aircraft response variable. 


To predict the mean rate N+(y) of exceedance crossings with 
positive slope of the aircraft response past the "level" y, one 

must compute the parameters , a* , ; therefore, as 

yz y z y s y s 

is evident from Eqs. (4,8) and (4.22), one must have available 
the power spectra $z(f) and the turbulence components 

o 

z(t) and Wg(t) in the model of Eq. (2.3). See Eqs. (4.17), 

(4.18), (4.20), and (4.21). These same two spectra are also 
required to predict the aircraft response probability density 
function p(y). Furthermore, computation of N+(y), using Eq. 

(4.8) or p(y) using Eq, (4,42), requires the probability density 
function p(bf) of the time-varying variance af(t) in Eq . (2.3). 

On the other hand, using Eqs. (4.28) and (4.46) to compute N+(y) 

and p(y), respectively, requires the central moments y^^.^of al 

O ^ X 

defined by Eq. (4.29). In addition, in deriving our expressions 
for N+(y) and p(y), it was assumed that the turbulence component 
Wg(t) in Eq . (2.3) was Gaussian. To check the consistency of 
this assumption, we require the probability density function 
p(ws) of Wg(t). A simple check of the Gaussian character of 
p(ws) can be made using a Gram-Charlier series representation 
of p(Wg), which requires the first few moments of ws(t). Finally, 
it can be shown that when the locally stationary response con- 
ditions of Eqs. (3-43) and (3.46) or Eqs, (5-18) and (5-22) are 
not satisfied, the power spectral density <I>g 2 (f) of the 

component a|>(t) is required to predict response exceedance rates 

and probability density functions. In summary, to predict the 

aircraft response statistics described above, one requires 

methods for computation of the power spectral densities 

(f), and $ 2 (f) of the turbulence components z(t), Wg(t), 
s ^ f 

and a^(t), respectively; methods also are required for computa- 
tion of the central moments and probability density functions 

y^^^and y4^\ and p(cff) and p(w ) of the turbulence components 


0 ^(t) and Wg(t). [The component z(t) in the model of Eq. (2.3) 
is, by assumption, Gaussian.] 
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In this section, we describe methods for computation of the 
above turbulence metrics from recorded turbulence velocities. 


Estimation of Power Spectra of z(t) and Wg(t) 

Let us first obtain an expression for the autocorrelation 
function R^^(t) using the model of Eq . (2.3). Multiplying w(t) 
by w(t+x) and taking the expected value of the resulting ex- 
pression yields 


E{w(t) w(t+r)} = E{w^(t) w„(t+x)} + E{w^(t) a™(t+T) z(t+T)} 

So S X 

+ E{Wg(t+T) a^(t) z(t)} + E{a^(t) a^(t+T) z(t) z(t+x)} 

= E{Wg(t) Wg(t+x)} + E{a^(t) a^(t+x)} E{z(t) z(t+x)}, (6.1) 


or 


R (x) = R (x) + R (x) R (x) 
w w a„ z 

s f 


( 6 . 2 ) 


where, in going to the second equality in Eq . (6.1), we have used 
the assumption that ws(t), Of(t), and z(t) are statistically 
independent and that z(t) has an expected value of zero. 


To interpret Eq . (6.2), let us examine Fig. 4, which was 
first discussed in Sec. 2. First, we note that the nominal 
correlation interval associated with the process af-(t) is much 
larger than that associated with z(t). To see this, examine, 
for example, the clumps of turbulence between approximately 
9 min 45 sec and 10 min 45 sec in Fig. 4. During this interval, 
the term Cf(t) in our model of Eq . (2.3) monotonlcally grows 
to a maximum from a small value and then decays monotonlcally 
back to a small value. The nominal correlation interval associated 
with Cf(t) is of the order of 1/2 min. On the other hand, the 
nominal correlation interval of the high-frequency process z(t) 
is of. the order of 1 sec. Thus, there is an order of magnitude 
or more separating the correlation scales of Of(t) and z(t). 

Since z(t) has, by assumption, zero mean value, this differenae 
between aorreZation scales means that Rqj,(t) is very nearly 


equal to 

negligible . 
may replace 


( 0 ) over the 
Hence, to a 


Eq. (6.2) by 


range of values of x where F^(xj is 
first and quite good approximation. 


not 

we 
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( 6 . 3 ) 


R„(t) - R„ (t) + R„ (0) R^(t) 

S I 

= (t) + aJ.-R^(T) , ■ 

S : 

where ai 5 R (0) Is the expecte(i value of the stationary process 

a|-(t). Let us now compare the nominal correlation intervals of 
the processes Wg(t) and z(t) by further- examination of Pig. 

It is evident that the correlation Interval associated with the 
low-frequency ( large-amplitude )■ process Wg(t) is at least of the 
same order of magnitude as that of the process af(t). Consequently, 
we may expect (i) to be very nearly constant and equal to 

(0) over the range of-Values of x where P ('xj is not negligible, 
s z 

These obser-vatlons are borne out by Fig. 5 , which is the 
autocorrelation function of the vertical, record shown in Fig. 4.* 
Equation (6.3) and the above discussion indicate that it-, .should 
be possible to approximate the autocorrelation function shown 
in Pig. 5 by a linear combination of R;^ (x) and R^-Ct) and that 

there should be an order “of magnitude or more difference between 
the nominal correlation Intervals associated with these two 
additive components. ,Thls situation is precisely what Pig.' 5 
Illustrates. In^fact, lf:we conceive ,of representing Rw_(t,) for 
X ^ 0 by a Maclaurln serlds, 

r; (o+) 

’’w ' '’w K . (6.1|) 

s s s 

then it is evident from Fig. 5 that the linear approximation 

gives a good representation of R;^_(x) out to a lag of about 

s 

10,000 ft (30^8 m) , whereas the quadratic approximation in Eq . 

(6.^) would appear to give a good approximation over the entire 
lag interval shown in Fig. 5. For lag values in the range of 
1000 to 10,000 ft (30^.8 m to 30^8 m), the deviations from the 
straight line shown in Fig. 5 may be attributed to statistical 


*Plgure 5 is the autocorrelation function expressed in terms of 
spatial lag 5 = Vx rather than temporal lag x, where V is the 
speed of the measurement aircraft. -This linear - scale factor 
on the ordinate has no effect on the concepts being discussed. 
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1000 ft 
(304.8m) 


VERTICAL GUST VELOCITY 
AUTOCORRELATION 


.LINEAR APPROXIMATION TOR^ (r) 
RwJt)«RwJO)+R'^ (0)|t| ® 


0 1 2 3 4 5 6 7 8 9 10 II 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 xIO^ 


5 (meters) 


FIG. 5. AUTOCORRELATION FUNCTION OF VERTICAL RECORD SHOWN IN FIG. 4. [MOUNTAIN 
WAVE CONDITIONS. AIRCRAFT SPEED 197 m/sec (646 ft/sec).] 







fluctuations In the estimate of R„(t) shown in Pig. 5 that are 
a result of the finite duration or the sample [which contains 
about 30 correlation intervals of the low-frequency component 

Ws(t) ]. 

In the lag interval of Fig. 5 from 0 to 1000 ft (0 to 30^.8 

m), we see superimposed on Rwg(T^)» as is predicted by 

Eq. ( 6 . 3 ). Examination of Eq. (6.3) at t = 0 indicates that, by 
setting R(0) = for both processes WgCx) and z(t), we have 


a 


2 

w 



+ 


a 


2 

z 




(6.5) 


where we have written 



( 6 . 6 ) 


which is consistent with the notation of Eq . (2.3). Examining 
the Intersection of the straight-line approximation with the 
ordinate in Fig. 5j we see that 



0.78 



0.22 


(6.7a,b) 


i.e., about 78 % of the contribution to the mean square value of 
w(t) in the vertical record is provided by the slow (low- 
frequency) component, whereas about 22% of the contribution is 
provided by the fast (high-frequency) component. These esti- 
mates are entirely consistent with the visual appearance of 
the vertical record in Fig. 4. 


The above comments immediately suggest a first approximation 
to R 2 (t). Solving Eq. ( 6 . 3 ) for R 2 (t), we obtain 

4 


s s 


( 6 . 8 ) 
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where. In the second line, we have substituted the Maclaurln 

series approximations to (t) given by Eq . (6,4) and where 

s 

the coefficients Rwg(0)> obtained either 

"by eye", as in the straight-line approximation of Fig. 5j or 
by a more sophisticated technique — e.g., least squares. 

Once Rj-Ct) Is solved, using Eq . (6.8), we may find $-(t) by 

forming the Fourier transform of R 2 (t). then be 

— O 

found by subtracting o|. From the power spectrum ^*yj(f) of 

the original record, as the Fourier transform of Eq. (6.3) would 
Indicate : 

s 

where a|. = R^(0) - since we have R 2 ( 0 ) = 1 by definition. 


Two minor refinements of the above procedure will now be 
discussed. First, we note that, in the high-frequency region 
where the contribution of O^^^(f) to 'J>;(f(f) is negligible, measured 

spectra consistently display a slope of -5/3, as predicted by 

the von Karman spectral forms. [The spectra (f) contaminate 

s 

iI> 2 ;(F) in the low-frequency region.] This observation suggests 
that we assume that 'J>z(F) the appropriate (transverse or 

longitudinal) von Karman form. Incorporation of the von Karman 
spectral form assumption for $z(f) eliminates the problems 
associated with statistical fluctuations in estimation of the 
functional form of 3>2(f). Estimation of ^>z ( f ) is then reduced 
to estimation of the Integral scale L of the component z(t), 
since a| = 1 by definition. 

A rough check of the von Karman form assumption for $ 2 (f) 
was performed as follows. Here, again, we have dealt with auto- 
correlation functions rather than spectra. The autocorrelation 
function dealt with is that computed from the vertic'al component 
of a sample case obtained from the NASA Langley MAT project, 
which is described in Ref. l8 . This autocorrelation function 
is displayed in Fig. 6 on the same scale as that shown in 
Fig. 5* However, it is evident from Fig. 6 that the relative 
contribution of the low-frequency component w„(t) is much less 
in this case than it was in the case displayed in Fig. 5. The 
time history from which Fig. 6 was computed is displayed as the 
top trace in Fig. 7- Comparison of the top traces in Figs. 7 
and 4 clearly shows that the relative contribution of Wg(t) in 
the vertical trace in Fig. 7 is much less than that in the 
corresponding trace in Fig. 4. 
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FIG. 6. 


1 


AUTOCORRELATION FUNCTION 
CONDITIONS. AIRCRAFT SPEE 










The autocorrelation function shown In Pig. 6 Is again dis- 
played on an enlarged scale In Fig. 8. To perform a rough test 
of the von Karman assumption using this autocorrelation function, 
we estimated (crudely) that ^ 0.10 for the autocorrela- 

tion function shown In Pigs. 6 and 8. Furthermore, for this 

rough test, we retained only the constant term Rwa(O) In Eq. 

s 

(6.4) as our approximation to R„g(r). Using this approximation, 
Eq. (6.8) reduces our estimate of Rz(t) to 


R,(t) = iR„(T) - ^R„ (0) 


( 6 . 10 ) 


When normalized, this approximation to Rz(t) is given by the 
continuous curve shown in Pig. 8 when considered as a function 
of the relabeled ordinate R 2 ;(t)/R 2 ;( 0 ). 


The encircled points shown In Fig. 8 are points of the 
(transverse) von Karman autocorrelation function plotted as a 
function of Rz(t)/R 2(0) for an Integral scale of L = 170.7 m 
(560 ft). We refer the reader to p. 253 of Ref. 19 for the 
mathematical form of the von Karman transverse autocorrelation 
function, as well as for a graph of this function. The encircled 
points shown in Pig. 8 were obtained by careful reading of points 
off the von Karman transverse autocorrelation function plotted 
on p. 253 of Ref. 19, after the points had been scaled to the 
Integral scale shown in Pig. 8. The fit of the circled points 
to the continuous curve Is easily within our reading error of 
the original curve in Ref. 19. Prom Fig. 8, we must conclude 
that the von Karman curve gives an excellent fit to the (small t) 
region of the curve that is not highly contaminated by the auto- 
correlation function Rw„(t) of the low-frequency component Ws(t) 

O 

in thie model of Eq . (2.3). 

A careful comparison of Figs. 5 and 6 Indicates that no 
measurable difference in the components Rz(t) in those two curves 
can be ascertained except for amplitude and integral-scale scale 
factors. Thus, it is likely that the assumption that Rz(t) and 
$z(f) have the appropriate von Karman forms is a good one. 

One further comment about Pig. 8 is in order. The von 
Karman form falls almost exactly on the original curve where it 
first touches the new abscissa. Thus, one could estimate the 
Integral scale of the component z(t) quite accurately in this 
case by using the point where it first touches the new abscissa 
(l.e., the lowest circled point). However, the Integral scale 
of the component z(t) in the model of Eq. (2.3) would be con- 
siderably overestimated if one computed it from the first 
crossing of the original abscissa (as is commonly done). 
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f (meters) 


FIG. 0. AUTOCORRELATION FUNCTION OF FIG. 6 ON EXPANDED SCALES. 
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The second refinement to the method for estimation of 
RgCx) given by Eq. (6.8) has to do with the assumption that 

Rap(T) = Rq^(O) = in the region of the t axis where Rz(t) is 

not negligible — l.e., the assumption that took us from Eq . 

(6.2) to Eq. (6.3). This assumption can be circumvented, because 
the normalized autocorrelation function 


R„ (T) 

A f 




( 6 . 11 ) 


of of(t) can be estimated using the turbulence model of Eq . (2.3), 
as we shall show shortly. Before showing how Po^(t) can be com- 
puted, let us first summarize the above described procedure for 
estimating $z(f) and $^^g(f), where we shall assume that Pcj^(t) 
is known. 


Procedure for estimating ^^(f) and Figure 9 il- 

lustrates the estimation procedure. The steps are: (1) Estimate 
a linear, quadratic, or possibly higher-order polynomial approxi- 
mation to Rwa^"^) neighborhood of t = 0, as is Illustrated 

o 

in Pig. 5. (2) Subtract this estimate of Rw„(t) from R„(t). 

The remaining function is our estimate of Ro-£.(t) R 2 (t), as is 

indicated by Eq. (6.2). Since, by definition, Rz(0) = 1, we 
have 


Of = 5' (6.12) 

i.e., we obtain an estimate of in this operation.* (3) Divide 
the estimate of Rq^(t) R 2 (x) by which we shall describe 

how to estimate shortly. This division yields an estimate of 
R 2 (t). (H) Determine the Integral scale that yields the best 

fit of the appropriate von Karman autocorrelation function to 
this estimate of R 2 (t). Note that Rz(0) will be unity. (5) Sub- 
tract the product of this von Karman estimate of Rz(x) and Op 
Pg (x) — i.e., Rap(x) R 2 (t) — from R^f(x). This subtraction 

yields a new estimate of Rwg(x). (6) Examine the behavior of 

*An independent alternative method for estimating cp = E(ap} is 
described in Appendix G. 
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this new estimate of (t) in the vicinity of t = 0 to ascertain 

if there is any obvious contribution from Raj«(T) Rz(t) remaining. 

If there is one, form a new estimate of (t) that eliminates 

s 

the contribution of Ro^(t) Rz(t) and repeat Steps (1) through 

( 6 ). (7) The final estimate of $z(f) is the appropriate von 

Karman spectral form with the integral scale provided by the 

above procedure. The final estimate of is the Fourier 

s 

transform of the final estimate of R;^g(T) determined in Step (5). 


Procedure for estimating Qa^(x), Here, it will be assumed 

that there exists a frequency f^ such that for all f ^ fo the 
contribution of the component Wg(t) of the model of Eq. (2.3) 
to the spectrum of w(t) is negligible in comparison with the 
contribution of the component wf(t) = Cf(t) z(t). The frequency 
ffl should be chosen on the low end of the frequency range of 
the portion of the spectrum that satisfies the - 5/3 decay 

law of the von Karman spectrum. This value of fo would cor- 
respond to an inverse wavelength of about 3 10 “^ cycles/ft 

(9.84 X 10-^ cycles/m) for the spectrum shown in Pig. 10, which 
is the wavenumber spectrum of the vertical record shown in Fig. 
4, whose autocorrelation function is shown in Pig. 5. 


The first step in estimating P 0 ^(x) from a turbulence 

record w(t) is to high-pass filter w(t), where the high-pass 
filter attenuates all frequency components of w(t) for f < f 
Denote the filtered record by Wh(t), as was done in Eq . (5.2 


Wj^(t) = A’ a^(t) Zj^(t) 


(6.13) 


where, here. A' is an arbitrary positive constant and z^(t) is a 
high-pass filtered version of z(t). Since Of(t) is nonnegative, 
by definition, we may express the absolute value of Wj^(t) as 

lw^(t)| = A' a^(t) |zj^(t)| . (6.14) 

Furthermore, since the record wj^(t) is available for manipulation, 
we may compute its absolute value and find the autocorrelation 
function of |wj^(t)|, which, according to Eq . (6.l4), is related 
to the autocorrelation functions of ap(t) and |z}^(t)| by 

= (AM“ , (6.15) 
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3.05x10'^ 3.05x10“® 3.05x10"* 3.05x10"^ 3.05x10"* 3.05x10"^ 

INVERSE WAVELENGTH, 1/X, cy/m 

FIG. 10. POWER SPECTRUM OF VERTICAL RECORD SHOWN IN FIG. 4. (CORRESPONDING AUTO- 
CORRELATION FUNCTION IS SHOWN IN FIG. 5.) [MOUNTAIN WAVE CONDITIONS, 
AIRCRAFT SPEED 197 m/sec (646 ft/sec).] 


\J\ 



where. In arriving at Eq. (6.15), the assumption that Of(t) and 
zj^(t) are statistically independent was used. Solving Eq . (6.15) 
for Raf(T) yields 



(t ) 




( 6 . 16 ) 


furthermore, dividing Eq. ( 6 .l 6 ) by R 0 „(O) — and using the defi- 
nition of Eq. (6.11) — yields 


(t) = 




{A')^ R (0) R| , (t) 
^f rh 


(6.17) 


It was shown in Ref. 5 that only rarely can we not assume that 
the autocorrelation function of wj^(t) in Eq . (6.13) is well 
represented by 

R ^ (t) = (A' )^ E{ap R^ (t) 

’"h ^ S 


= (A»)" R_ ( 0 ) R (t) . ( 6 . 18 ) 

®f ^h 


Equations (6.17) and ( 6 . 18 ) provide us with the means to estimate 
the normalized autocorrelation function P 0 ^(r). The numerator 

of Eq . (6.17) can be computed directly from the absolute value 
of the filtered turbulence record. Also, we can compute Rw^(t) 

from the same (unrectlf led ) record. Furthermore, since z(t), 
by assumption, is a stationary Gaussian process, zj-i(t) also is 
stationary and Gaussian. We therefore can compute the quantity 
(A’)^ Raf(O) from the measured autocorrelation function 

Rwh(T) by using the well known relationship [e.g. , p. l 66 of 

Ref. 8 ] 

2 I I 

(A<)‘ R (0) H (T)' . ( 0 )]> . [R (T)]“ 

f I h I ' h h 

+ R (t) arccos C-R„ (t)/R ( 0 )] - 5 R„ (t)| . (6.19) 

% ''^h % '^h ) 
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Thus ^ aaaording to Eq . (6.17) ^ Pojy('^) ^ccy be estimated by 

dividing the autoaorvetation function R^u-^^(x) \_whioh is to 

be computed directly from the absolute value \w-j^(t)\ of the 
filtered record W}i(t)~\ by the right-hand side of Eq. (6.29) 
[which is computed from Tt ,) ] . The autocorrelation function 

computed directly from the (unreatified) record 

w-^(t) .* 


We may check the 
Ing Eq. (6.17) at t = 
R|Whl(0) = E{|wjJ^} = 


consistency of the above method by evaluat- 
0 and T = “>. First, consider t = 0. Since 
E{w^}, Eq. (6.17) yields, at t = 0, 



( 0 ) 


E(w^} 


(A’ )^ (0) R „ (0) 

'"f |"hl 


( 6 . 20 ) 


*The approximation made In Eq . (6.l8) can be avoided If the 
"arcsln law" Is employed to compute R|Zj^|(t) In the denominator 

of Eq . (6.17). We used the arcsln law successfully with turbu- 
lence data In work reported In Ref. 6 (e.g., pp . 16-21). In 
using the arcsln law In the computation of Pq^(t), one must 

replace the denominator. of the right-hand side of Eq . (6.17) 
by 


f R^ (0) ! (t) arccos[-R^, (x)] - R^ , (x)} , 

h ( T h h h h ) 

where 

R„, (x ) = sln[^ Rp (x ) ] , 

^h ^ 

where Rq(x) is the autocorrelation function of the hard-clipped 
version of.wj^(x) that Is formed by setting the hard-clipped 
waveform equal to +1 where it is positive and -1 where It is 
negative . 
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The denominator of Eq. (6.20) must be evaluated from Eq. (6.19); 
for T = 0, using the fact that arccos(-l) = it, we find from 
Eq. (6.19) that 

(A')= R„^(0) - 1} ” 5 

= (0) = E{w^} ; (6.21) 


hence, by combining Eqs. (6.20) and (6.21), we have Pq„( 0) = 1, 
which is, of course, correct. ^ 

Now let us evaluate Eq. (6.17) at t = “. Recognizing that 
R|wh|^“) = {E[|wj^|]}^, we see that Eq. (6.17) yields, at t = <», 

{E[|w, |]}2 

P (oo) = , (6.22) 

(A*)"" R„ (0) R ^ (oo) 

\^h\ 


Evaluating the denominator using Eq. (6.19) and recognizing that 
Rw^(“) = O, we have 




Combining Eqs. (6.22) and (6.23) yields 

{E[ |w^| ]}^ 


(“) = ^ 


IT h' 


However, from Eq . (6.1^4), we have 


E[|w^|] = A’ E[a^] E[|z^I] , 


whereas from Eq, 


E[w"] 


(6.13) we have 
(AM" E[aJ] E[z^] 


(6.23) 


(6.2^4) 


(6.25) 


( 6 . 26 ) 
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Substituting Eqs. (6.25) and (6.26) into Eq. (6.24) yields 


{ E [ a ^]}2 { E [| z ^ 1]}2 


(6.27) 


However, Zh(t) is a Gaussian process with zero mean value; hence, 
we have [e.g., Eq. (4.4-25) on p. l66 of Ref. 8] 


E[ 



Therefore, Eqs. 




(6.27) and ( 6 . 28 ) yield 
{ E [ a ^]}2 

E[ap 


( 6 . 28 ) 


(6.29) 


which is the correct value. The method of Eqs. (6.17) and (6.19) 
therefore checks, exactly, in the limiting cases t = 0 and t = «>. 


Estimation of Power Spectrum of a|(t) 


A procedure similar to that described above can be used to 
estimate the autocorrelation function Ral(T) of Of(t). Squaring 
Eq . (6.13), we have 


w^(t) = (A')^ o^(t) z^(t) . (6.30) 

From the assumed statistical independence of ap(t) and z^(t), 
it follows that 


R„2(t) = (A')** R„2(t) R„2(t) 
% ^h 


(6.31) 


where we have defined 



(t ) 


E{w^(t) w^(t+x)} 


(6.32) 
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with comparable definitions holding for R 0 i(x) and Rz^Ct)- 
Solving Eq. (6.31) for Rq^Ct) yields 



(t) 


R^2(t) 

^h 

(A')" R^2(t) 
^h 


(6.33) 


By squaring a given high-pass filtered recording wjj(t), we 
can compute R^ 2 (x) directly. To compute Rz^(''^)j we begin with 

the approximation of Eq. (6.18), which. can be expressed as 

R„ (t) = (A')2 E[c2] r (t) , (6.34) 

^h ^ ^h 

which is equivalent to the expression 


E{w^(t) w^(t+x)} :: E{[A'/E(ap Zj^(t)] [A’/E(a|.) z^(t + x)]} .(6.35) 

Now, if Eq. (6.35) were an exact expression^ it would follow that 

R„ 2 (x) = (A')"* {E[a"]}2 R^ 2 (x) , (6.36) 

% ^ ^h 


from which it would follow that 


R ..2 (x) 


R , ,2 (x) 

^h 

(A*)** {E[ap}2 


(6.37) 


But, from Eq . (4.4-26) on p. l 66 of Ref. 8 , we have for stationary 
Gaussian processes {w^(t)}. 


R^ 2 (x) = [R^ (0)]" + 2[R^ (X)]" . (6.38) 

Substitution of Eq. ( 6 . 38 ) into Eq. (6.37) yields 
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9 


(6.39) 


R„2(t) 


CR„ (0)]2 + 2[R (t)]2 

% 

(A')“ {E[ap}2 


the validity of which depends on the stationary Gaussian assump- 
tion for (zh(t)} and' the approximation of Eq. (6.34). Combining 
Eqs. (6.33) and (6.39) yields, finally. 


R^2(t) « 


{E[a^]}' 




[ R „ ( 0)]2 + 2 [ R ^ ( t )]= 

^h "h 


(6.40) 


which is the desired expression.* Notice that R^j^(t) can be 

*The approximation of Eq. (6.34) can be eliminated in deriving an 
expression for Rrr 2 (T) using the "arcsin law" as was the case with 

f 

P(jf.(T) considered in the previous footnote. Prom Eq . (6.33) > we 
have, exactly, 

R^2(t) R„2(t) R^2(0) 

o /-N A ^f _ % ^h 

R 2(0) RTToT R^p(t) 

f Wj^ 

Applying Eq . (6.38) to the Gaussian process z^(t) yields 
R^2(t) = [R^ (0)]"' {1 + 2[p^ (t)]M . 

Using the "arcsin law", we have 


p (t ) = sin[^ Rq (t ) ] , 

^h 

where Ro('c) was defined in the last footnote. The above ex- 
pressions can be evaluated to yield p„ 2 (t) and, therefore, 

Of 

R„ 2 (t) in terms of Rrf2(o). To compute R„2(0), we note that 

O £» ^ f* 

R^2(“) {E[ap>2 

which can be solved for R^j 2 ( 0 ) in terms of p^j 2 (<») and {E[a^]}^. 
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evaluated directly from the high-pass filtered turbulence . 
sample Wj^(t), whereas can be evaluated directly from 

the squared sample w^(t). E{a|.} can be evaluated by the method 


Illustrated in Fig. 9 or by the method described in Appendix G. 
Notice also that, from the definition of R„ 2 (t), we must have 

vJ ^ 


lim 


R 2(t) 

Of 


{E[op}" 


(6.41) 


To check the consistency of Eq. (6.40), we shall evaluate 
it for the limiting cases t = 0 and t = “. For t = 0, Eq. (6.4o) 
yields 


R„2(0) = 

Of 


{E[op}2 




{E[w2]}2 


2{E[wJ]}2 


= { E [ o 2]}2 


E{w^} 

3{E[w^]}2 


(6.42) 


However, from Eq . (6.13), using the statistical independence of 
o^(t) and z^(t), we have 


E{wJ} = (A’)** E[oJ.] E[z^] (6.43) 

and 

E[wJ] = (am" E[o^] E[z2] . (6.44) 


Substituting these quantities into Eq. (6.42) and cancelling 
identical quantities in the numerator and denominator yield 

E[z;] 

R^2(0) = E[oM . (6.45) 

""f ^ 3{E[z^]}2 

However, zj^(t) is assumed to be a stationary Gaussian process 
with zero mean value. It follows (e.g., p. 221 of Ref. 10) that 
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5 


(6.46) 


E[zJ^] = 3 {[zJ]> 


hence, we have from Eqs. (6.45) and (6.46) 

R„2(0) « Elo'l-] , 

°r ^ 


(6.47) 


which, of course, is an Identity. Thus, Eq. (6.40) yields the 
correct value at t = 0. 


For T = “, it follows directly from Eq. (6.40) that we may 

write 


R 2 (~) 

Of. 


{E[ap}2 


{E[w^]>" 

{E[w^]}" 


= {E[ap}'‘ 


(6.48) 


which is the correct limiting value given by Eq. (6.4l). Con- 
sequently, Eq. (6.40) provides the correct exact limiting values 
at T = 0 and t = «>. 


To obtain the power spectral, density ^* 52 (f) of Of(t), we 

Fourier transform Eq. (6.40). Recognizing Eq. (6.41)j we may 

express <I'^ 2 (f) as 
Of 


$„2(f) 

°f 


= {E[a|]> } 6(f) 


R,.,2(t)-[R, (0)]^- 2[R, (t)] 


w 


h 


w 


w. 


-«>\ [R„ (0)]" + 2[R„ (t)] 


w 


h 


w. 


-i27TfT V 
e dT>, 


(6.49) 


where 6(f) is the Dirac delta function and where may 

be evaluated using the method described in Appendix G or from 
the method illustrated in Fig. 9. All that are required to 
evaluate the integrand ir^ Eq. (6.49) are the autocorrelation 
functions of wji(t) and w^(t)^ both of which can be evaluated 
directly from the high-pass filtered turbulence sample wh(t) . 

The Fourier transform in Eq. (6.49) would, of course, be carried 
out numerically . 


83 



Estimation of Moments and Probability Density of Of(t) 


Moments of Of. Let us now turn to the development of 
methods for determination of the moments and probability density 

function of o%(t) . The moments yfa^required for computation of 

X ^ f* 

the response exceedance rates N+(y; and probability density 
function p(y) In Eqs . (4.28) and (4,46)^are the central moments - 
l.e., moments taken about the mean of Of. See Eq. (4.29)* 

Here, we shall develop a method for determining the moments 

taken about the origin — l.e., 

^f 


M 


(k) A 
^2 - 


.00 

(op^ p(cpda^ 


0 


E E{[ap^} 


(6.50a) 


(6.50b) 


It is easy to show [e.g., Eqs. (98) and (99) on pp . 273 and 274 
of Ref. 16] that the nth central moment can be computed from the 

moments k = 0,l,2,...,n by the formula 

Of 


(n) 
a 



, n = 1,2,3»*»* > 


(6.51) 


where we note that 
function; hence. 


( 0 ) 
MS ^ 

Of 


Is the 


area under the probability density 



(6.52) 


and 


n \ _ n! 
k I ” k ! (n-k ) ! 


are the binomial coefficients. Notice, 
n = 1 and 2, Eq . (6.51) yields 


V 




(6.53) 


in particular, that for 
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(6.54a,b) 



Equation (6.5^a) states that the first central moment Is zero, 
whereas Eq. (6.54b) Is the familiar expression for the second 
central moment, which Is the variance. [See Eq. (4.23).] 

Let us turn now to developing a method for determining the 

moments n = 1,2,3,... defined by Eq. (6.50). To do so, 

we shall again assume that we have available the high-pass 
filtered turbulence sample w^(t) described earlier, which has 
the form of Eq. (6.13), where A’ is an arbitrary positive constant 
and z^(t) Is a high-pass filtered version of the Gaussian turbu- 
lence component z(t) defined by Eq. (2.3). Let us now write an 
expression for the nth moment of w^(t), which Is expressed In 
terms of a|.(t) and z^(t) by Eq. (6.30). Using the fact that 
Of(t) and z(t) are assumed to be statistically Independent, It 
follows immediately that we can write the nth moment of Eq. 

( 6 . 30 ) as 




from which we 


may solve for E{[a^]^} as 


E{[a^]’^} 


E{[wJ]''} 


(6.55) 


( 6 . 56 ) 


Since we have assumed that the high-pass filtered sample w;^(t) 
Is available, we can numerically square it and then compute the 
first several moments E{Cw^]^^}, n = 1,2,3,-... Furthermore, 
since z(t) is assumed to be Gaussian, Z)-^(t) also is Gaussian; 


hence, we can compute the moments E{[z^]*^} in terms of E{z^}, 
as we shall now show. 


It Is known (e.g., pp . 233-236 of Ref. 9) that the pro- 
bability density function of the square of the Gaussian variate 
z^ Is given by the chi-square density with one degree-of- freedom. 
Letting 



and 

= E{z^} 


(6.57) 

( 6 . 58 ) 
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we have, from p. 236 of Ref. 9, for the probability density of 


P..2(X) 

“h 


1 


|/2ir a 


^-x/(2aM 


X > 0 


X < 0 


(6.59) 


To determine an expression for the moments of Zhj 


M 


^2^ E E{Czn'^} = 


h- 


M 


(n) 


( 6 . 60 ) 


we shall want to consider the normalized variable, which Is de- 
fined from Eqs. (6.57) and ( 6 . 58 ) as 





( 6 . 61 ) 


Let P^(0 denote the probability density of C. 


Then, since 


p 2 (x)dx = p^ (C )d^ 


( 6 . 62 ) 


It follows Immediately that the moments 
from the moments of C by 


of z^ can be computed 




p 2 ( X ) dx 
^h 


.00 

)” p^ (5 )d^ 

0 


( 6 . 63 a) 


= , ( 6 . 63 b) 

( n ) 

where we have used Eqs. ( 6 . 6 I) and (6.62) and where denotes 

the nth moment of the variable C = x/a^. 

Using Eqs. (6,.6l) and (6.62), it follows from Eq . (6.59) 
that the probability density of ^ is 
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1 




£ > 0 




=< 

(o , 5 < 0 , (6.6H) 

where we have used the fact that (dx/dC) = . Equation (6.64) 

is the usual normalized form of the chi-square density function 
with one degree-of- freedom. The moments of 5 are given on p. 405 
of Ref. 20. See Eqs. (65) and (66) of Ref. 20, using m = 1. 

These moments are 



r(|+n) 


,n 


(6.65) 


where we have used the fact that F(l/2) = /it. Using the relation- 
ship r(n+l) = n r(n), together with r(l/2) = /F, one may 

successively write out F(^n) for n = 1,2,... to discover that 

F(^n) = ^ ^ ^ ^ . (6.66) 

Combining Eqs. (6.63b), (6.65), and (6.66) yields the desired 
expression for : 


= 1-3-5 ... (2n-l)a^" , (6.6?) 

A Zp, 


where we have Inserted the subscript Zj, on the a to denote that 
it represents the standard deviation or the variable z^. Com- 
bining Eqs. ( 6 . 56 ) and (6.67) and using the notation ' of . Eq . ( 6 . 6 O) 
together with the definition 


E[w^] = ol , (6.68) 


we obtain for the nth moment of a|. 
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( 2 n-l) 


9 


E{[a2]’^} 




(A'a, 




, 2 n 


1-3-5 


n = 1 , 2 , 3 ,... . 

(6.69) 


To put the above expression in a more useful form, we note that 
for n = 1, Eq. (6.69) becomes 


E[wn 

E[ 0 ^] = . (6.70) 

(A'a^ 

^h 


Solving Eq. (6.70) for (A'Ozvj)^ and inserting that expression 
into Eq. (6.69) yields, finally. 


E{[aJ]’^} = M 



j ECaJ.])'" E{[w2]>^} 

( EEwJ]) 1 . 3.5 ... ( 2 n-l) 


n 


1 2 

i 3 • . . 


(6.71) 


Equation (6,71) is the desired expression for the moments 

mLzj n = IjBfSj,.,, of Equation (6.51) yields the ae-ntral 

moments v^Lijrom these moments. In evaluating the right-hand 
f 

^ J 2 >2 % 

side of Eq. (6.71) y the moments }, n = 2,2,3,... are to 

be computed directly from the squared high-pass filtered turbi^- 
lence sample Wf^(t) described by Eq . (6.1Z). The quantity 

can be determined from the method illustrated in Fig. 9 or the 
method described in Appendix G. 

Probability density function of o^. :^et us now turn to 
estimation of the probability density of Cf. Since we already 
have an expression, Eq. (6.71), for the moments of Of, it will 
be convenient to develop an approximation for the probability 
density function of Of from these moments. To do so, we shall 
use an extension of the Gram-Charller series. 


The Gram-Charller series is an expansion of a probability 
density function, the first term in the expansion being a normal 
probability density function with the correct mean and variance 
and the coefficients of the, say, N correction terms being 
chosen so that the moments through order N of the series approxi- 
mation are set equal to the moments of the original distribution. 
The Gram-Charller series is described in this manner on pp . 270- 
272 of Ref. 14. However, the Gram-Charller series usually is 
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not derived using this moment equivalence criterion — e.g.. 

Ref. 9, pp. 222-223 or Ref. 20, pp . 136-137. Using the moment 
criterion for choice of the correction term coefficients, this 
writer has extended the Gram-Charller series to completely 
arbitrary "base density functions" (analogous to the normal 
probability density) on pp. 269-278 of Ref. l6 (where the 
application was to a different problem). 

Since the random variable under present discussion cannot 
be negative, the normal probability density function is not an 
appropriate "base density function" for expansion of a|>. However, 
the gamma density function — e.g., pp. 220-221 of Ref. 10 — is 
appropriate since it has the flexibility required in the present 
application, yet does not permit negative values of the variate 
It describes. The gamma probability density has the form 

(r^ e"^^ , V > 0 

p(v;y,?^) =j 

(o , V < 0 , ( 6 . 72 ) 

where F(y) Is the gamma function and where V = a|. is the variate 
being described. 

Notice that p(V;y»^) contains two "free" parameters, y and X. 
The mean and variance of p(V;y,X) are related to y and X by 


E{V} = M. 


(1) - Y 


V 


E[V^] - {E[V]}^ = M. 


X 

,( 2 ) 

V 


^ x^ 


(6.73) 

(6.74) 


from which we may solve for y and X In terms of the mean and 
variance to give 


Y = 


X = 






(6.75) 


(6.76) 
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Of particular Interest, In the present application. Is the 
capability of p(V;y,A) to approximate the density function of 

Y = Of when the relative standard deviation (or relative 
variance) of Of Is small — l.e., when the ratio of the standard 

deviation of a ^ to the mean of is' small. According to Eq. 

(6.75), this ratio (for the variate V) is given by Thus, 

the asymptotic form of p(V;y»^) for large y is of interest. 
However, it is immediately evident from the form of the charac- 
teristic function of the gamma density p(V;y,?^) — e.g., p. 221 
of Ref. 10 — that, according to the central limit theorem, the 
gamma density function approaches a normal or Gaussian density 
as Y furthermore, in the actual limit y = °°j p(V;Yj^) 

becomes a Dirac delta function located at V = E{V} = when 

Y and X are chosen to satisfy Eqs. (6.75) to (6.76). Thus, for 
small values of relative standard deviation, p(V;y,^) has the 
general appearance of a Gaussian density centered at E{V>, but 
with a truncated tail so that the density is zero for negative 
values of V. For large values of relative standard deviation, 
p(V;Yj 1) can take on a considerable range of shapes, as an 
examination of Eq. (6.72) will indicate. Furthermore, it is of 
some Interest to note that, if Of were to be exactly normally 
distributed with zero mean value (which we have ruled out by 
hypothesis), then the probability density of Of would be exactly 
described by Eq. (6.72) with y = 1/2. 


From the above comments, we might expect Eq. (6.72) to 
provide a reasonable approximation to the probability density of 
V = of* if the mean and variance of the density are chosen to 
take on the mean and variance of the variate V = a|., a decision 
which would require that y and X be chosen by Eqs. (6.75) and 

( 6 . 76 ). Thus, computation of only two moments and 

o ^ ^ f* 

by Eq . (6. 71) should provide a reasonable first estimate of the 
probability density function of a^. 

Let us now consider the extension of the Gram-Charller 
series to the base density function given by Eq . (6.72). This 
extension is provided on pp . 272, 273j and 276 of Ref. I 6 , where 
in Eq. (II 8 ) of Ref. I 6 , we must set A = X/f(Y), since the 
probability density p(V;Yj^) must have unit area. It follows 
from Eqs. (97), (ll4), and (II 8 ) of Ref. I 6 that our generalized 
form of the probability density of V = can be written as 


rTTT ® 


N 

1 + y b’ p (V) 
n=3 


P(j2(V) = 

Of 


V > 0 


, V < 0, ( 6 . 77 ) 


90 



where 


= iFi(-n;Y;AV) ( 6 . 78 a) 

n (-n)i. V 

= I r" , n = 1,2,3,... ( 6 . 78 b) 

k=0 k!(y)j^ 


are proportional to the generalized Laguerre polynomials, where 
iPj(-n;y;XV) is the confluent hypergeometric function, and where, 
in Eq. ( 6 . 78 b), we have used the notation 


(1^)0 = 1 

(u)j^ = u(u+l) ... (u+k-1) , k > 1 . (6.79) 


The expansion coefficients b^ In Eq. (6.77) are related to those 
in Eq. (97) of Ref. I 6 by 


b’ 


n 


^ = r(Y) 

A X 


b 


n 


( 6 . 80 ) 


as may be seen by examination of Eqs. (97), (ll4), (II 8 ), and 
(122) of Ref. 16 , and by recognition of the fact that, in the 

present application, we must have 1. Consequently, by 

substitution of Eqs. (9^b) and (117) of Ref. I 6 into Eq. (6.80) 
above, we find for our expansion coefficients b’. 


b’ 

n 


(y) 


n 


n! 


n 

I 

k=0 


a 


nk 


^f 


( 6 . 81 ) 


where, from Eq. (II 6 ) of Ref. I 6 , we have 




( 6 . 82 ) 


which are the coefficients of the polynomials defined by Eq. 
( 6 . 78 b) above. In Eqs. ( 6 . 8 I) and (6.82), we have used again 
the notation of Eq. (6.79). 
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Equations (6,77)^ (6.78)^ (6.81)^ and (6.82) aotleotiveZy 

describe the desired approximation to the probability density 
function of where the two parameters y X are to be 

evaluated by Eqs. (6.75) and (8.76) from the first two moments 
of V = Ojr. The series in Eq . (6.77). can be truncated at any 

value oj the index N. According to Eq. (6.81) ^ for any integer 
value of N ^ 2j we can evaluate the expansion coefficients b^ 


. (6 '.77) from the sequence of N moments 

that are to be determined by Eq- (6.71) 


Moreover i 


Eq. 


in Eq 

(6. ,77) has the property that we may add additional terms without 
changing the values of the coefficients of the terms previously 
determined. It has been shown in Ref. 16 that, for any value of 

(1) ,.(2) „(N) 

J 


N ^ 2^ the moments M y M ^ ...j Af of the approximation 

f(j^(V) given by the right-hand side of Eq. (6.77) are set equal 

expansion co- 


7 


to' the moments M^^^ , when the 


efficients are determined by Eq. (6.81) and when y and X are 
determined by Eqs. (6.75) and (6.76). 


For values of V large relative to the mean E[V] = y/X, the 
tall In the approximation to P 02 (V) given by Eq. (6.77) can go 

negative for N > 2. Some judgment will have to be used In 
choosing a value of N to prevent this occurrence and possibly 
other undesirable behavior of the series of Eq. (6.77). It Is 
unlikely that much accuracy of any utility will be gained by 
using values of N larger than ^ (two correction terms In the 
right-hand side of Eq. 6.77), and, for practical purposes, the 
approximation given by Eq . (6.72) (no correction terms In the 
right-hand side of Eq. (6.77) will probably be adequate In most 
cases . 


Estimation of Moments and Probability Density of Wg(t) 

Moments of Wg. To form our estimate of the probability 
density Pwg(wg) of Wg(t), we shall use the moments of Wg(t) to 

generate the coefficients of a Gram-Charller expansion of PwcC^s)* 

Moments will be used here because the moments of Wg(t) are 
particularly easy to compute from the moments of Wh(t) and w(t). 
For the present application, we shall require the moments of 
w^(t), rather than the moments of Its square w^(t). 
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According to the turbulence model of Eq. (2.3), the tur- 
bulence record w(t) is the sum of the "slow" and "fast" pro- 
cesses Ws(t) and Wf(t); l.e.. 


w(t) = Wg(t) + w^(t) 


(6.83) 


The moments 


M 


(n) _ 


= E[w^], n = 1,2,3, 


w 


(6.84) 


of the record w(t) can be computed directly, 
for now, that we have available the moments 


Let us assume. 


M. 


= E[wJ], n = 1,2,3,... 


(6.85) 


of the "fast" component Wf(t), It is shown in Appendix H that 
the moments of the "slow" component 


M 


= E[wg], n = 1,2,3,... 


( 6 . 86 ) 


can be computed sequentially from the sequences ^ and , 

n = 1,2,3,... by the relationship 


M 


(n) 


w. 


= - Jo (s) 


where 


jvi(O) = p^(0) = 2 

"^f ^s 


( 6 . 88 ) 


are the areas under the probability density functions of Wf(t) 

and Ws(t) and where (^) are the binomial coefficients defined by 

Eq. (6.53). The first four of the sequence of the relationships 
of Eq. (6.87) are 
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Wg w 


w^ w ^ w^ w^. w^ 


iyi^3) _ ]vi(3) _ []vi(3) ^ 3iy[^^^ 
w w w« w^ w w„ w 

S 1 I s I s 


M 


w w w„ w„ w w„ w w„ w„ 


f 


f 


(6.89) 

Notice from Eqs . ( 6 . 87 ) and ( 6 . 89 ) that to compute the moments 
through order N of Wg(t), only the moments through order M of 
w(t) and wf(t) are required. 

The moments n = 1,2, 3,-.. of the component wf(t) are 

to be computed from the high-pass filtered version w^Ct) of the 
turbulence record w(t), where this high-pass filtered waveform 
w^(t) was discussed earlier in connection with Eq . ( 6 . 13 ). It 

is shown in Appendix G that the moments n = l,2,3j... can 

/ \ ^ 

be computed from the moments , n = l,2,3>-«* of the record 

Wh(t) by h 


I^(n) ^ j^n j^(n)^ ^ 1,2,3,... 

Wf w^ 


(6.90) 


/ N 

where the moments are, of course, defined by the relationship 


i E[w[J] , n = 1,2,3,... . 

h 


( 6 . 91 ) 


The positive constant K in Eq. (6.90) is shown in Appendix G to 
be given by 

,00 


K = 


S^(f)df 


H^(f) i ^ S^(f)df 


( 6 . 92 ) 
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where 32 (f) Is the power spectral density of the component z(t) 
in the model of Eq. (2.3) and where Hh(f) is the high-pass 
filter complex frequency-response function. As we described 
in Sec. 6.1, we shall generally want to assume that 32 (f) is 
the appropriate von Karman spectral form. Moreover, since 
E[z^] = 1 [according to Eq. (2.3)]j the numerator in the right- 
hand side of Eq. (6.92) will be unity. 

Once the integral scale of z(t) has been determined and the 
high-pass filter has been chosen, the constant K can be computed. 
Furthermore, from a turbulence record w(t), we may generate the 

moments and the high-pass filtered record Wh(t), from 

which we may compute its moments Prom these moments and 

the constant K, we can use Eq. (6.90) to compute the moments 

which we may then combine with the moments using 

” I ^ r> ^ ^ 

Eqs. ( 6 . 87 ) or (6.89), to compute the moments M^ . 

s 


Probability density function of Wg. To generate an estimate 
of the probability density of Ws(t) using the Gram-Charlier 
expansion, we must use the central moments of Ws(t). If w(t) 
has zero mean value, then it is evident from the first of the 


four equations in ( 6 . 89 ) that we should have = 0; l.e., 

fl^ (n) ® 

when =0, the moments M„ , n = 1,2, 3, . . ., are the central 

S S 

moments. We shall assume that = 0 in the following dis- 

s 


cusslon . 


It is shown on pp . 270-272 of Ref. 14 that the Gram-Charlier 
expansion of a probability density function, say, PWg(Wg), may 
be expressed as 



where, here, 

a = {E[w|]}^ = (6.94) 

”s 

and where only two correction terms have been retained, the first 
being an odd function of Ws and the second being an even function 
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of Wg . These two correction terms are adequate to determine, 
how closely Pwc-^^s^ conforms to a Gaussian probability density 
function. It may be shown that the results provided on pp . 
277-278 of Ref. 16 also lead to the result given by Eq. ( 6 . 93 ). 

To evaluate the parameters in Eq. (6.93), one requires cr = 

, and These quantities are to be, evaluated using 

Eq. (6.89). 
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APPENDIX A 

DERIVATION OF REQUIREMENT FOR NEGLIGIBLY SMALL CORRELATION 
COEFFICIENT BETWEEN A NONSTATIONARY PROCESS AND ITS DERIVATIVE 

Consider an arbitrary function y(t) and Its derivative y'(t) 
Integrating the product of y and y’ by parts, we find 


_1_ 

At 


’t+At 

t+At 

rt+At 

y(?)y'(5)d? = ^ y"(?) 

- 

y(E) y'(?)d5 

t 

t 

t 


r 

At 


r 

•t+At 

y^(t+At) - y*(t) - 

y(E) y’COdE 

• 

t - 


(A.l) 


hence, we have, exactly, 
ft+At 


JL 

At 


y(5) y’(C)d? = 


1 y^(t+At) 

2 At 


- y"(t) 


(A. 2) 


Let us now assume that {y(t)} Is a generally nonstationary sto- 
chastic process. Taking the expected value of Eq. (A. 2) and then 
Interchanging the order of expectation and Integration operations 
on the left-hand side gives 


At 


ft + At 

E<y(5) y'(5)>dE - | . 

t . 


(A. 3) 


If we now take the limit At -► 0 In Eq. (A. 3), we have, assuming 
that E{y(?) y'(?)} Is continuous, 

E{y(t) y’(t)} = I ^ E{y*(t)} , (A.H) 

which holds for nonstationary and stationary processes. For 
wlde-sense stationary processes, the right-hand side of Eq. (A. 4) 
Is zero; hence, for wlde-sense stationary processes, Eq. (A. 4) 
reduces to the usual result E{y(t) y'(t)} = 0. 
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Denoting by by(t) and Oyt(t) the mean-square values of 
y(t) and y'(t) and assuming that 


E{y(t)} = 0 


for all t 


(A. 5) 


we may express the correlation coefficient of y(t) and y'(t) as 


n (f) = E{y(t ) y ' (t)} 

Oy(t ) Cy, (t) 


1 
2 


ay,(t) cjJ(t) 


«y 


(A. 6) 


which Is a completely general result for nonstationary processes. 

At this juncture, we shall assume that the locally stationary 
conditions of Eqs. (3.^1), (3.^3)» and (3.^6) are satisfied; 
hence, the response-process conditional instantaneous spectrum 
Is well approximated by Eq. (3.^0). When these conditions are 
met, we shall show that the right-hand side of Eq. (A. 6) Is 
negligibly small In comparison with unity. In carrying out the 
heuristic proof to follow, we shall further assume that the 
"slow" component of excitation Ws(t) In Eq. (2.3) Is zero, which 
implies that $^^(f) In Eq . (3.^0) Is zero. Removal of this 

stationary component of the response process has the effect of 
Increasing the magnitude of the right-hand side of Eq . (A. 6), 
which Is zero for stationary processes; hence, the assumption 
that $Wc,(i’) is zero is conservative. 

Integrating Eq . (3-^0) over -oo < f < oo and using a funda- 
mental property of the Instantaneous spectrum — e.g., Eq . (12a) 
of Ref. 7 or Eq. (2.4a) of Ref. 6 — we have 


(t ) 

y 


fCO 

$y (f ,t |a^)df 

t 

^CD 


»» a^(t)aj 
^z 


(A. 7) 
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where we have set equal to zero In Eq. (3.^0) before 

carrying out the above Integration, and where $v Is the mean- 

J z 

square aircraft response to the stationary component z(t) of 
the "fast" turbulence component 'oftt) z(t) In Eq. (2.3); l.e.. 



-CO 


$^(f) 


|H(f)|^df . 


(A. 8) 


Equation (A. 7) Is the desired expression for cXy(t) for use 
In Eq. (A. 6). Let us now obtain an expression for ffy»(t). The 
locally stationary response approximation of Eq. (3.^0) implies 
that we may write the response process as 


y(t) s a^(t) y^(t) 


(A. 9) 


where yz(t) Is the aircraft response to the component z(t) of 
Eq . (2.3). Differentiating Eq. (A. 9) gives 

y'(t) = a^(t) y^(t) + a^(t) y^^Ct) ; (A. 10) 

hence , 

a^,(t) ^ E{[y'(t)]M 

= a^(t) E{[y^]^} + 2a^(t)a^(t) E{y^y^} + [a^(t)]^E{yp 
= a^(t) o^, + [aL(t)]^ , (A. 11) 

^ rr ^ ^ r7 


where the middle term in the second line of the above equation 
is zero because E{y^yz) = 0, which follows from the fact that 
the process {yz(t)> Is stationary. Now, for locally stationary 
processes that satisfy Eq. (3*^0), we would expect the second 
term in Eq . (A. 11) to be small In comparison with the first 
term. Neglecting the second term yields the locally stationary 
approximation to a^,(t) — l.e., 

y 

aj,(t) = a^(t) aj, ; (A. 12) 

y ^ z 
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hence, from Eqs. (A. 7) and (A. 12), we have 


a (t) 
""y’Ct) 



(A. 13) 


which Is Independent of time. Moreover, we note that the approxi- 
mation given by the right-hand side of Eq. (A. 12) Is always less 
than the right-hand side of Eq. (A. 11); hence, the approximation 
of Eq. (A. 12) has the effect of Increasing the size of the right- 
hand side of Eq. (A. 6); l.e.. Insofar as the present discussion 
Is concerned, the approximation of Eq . (A. 12) Is conservative. 

We may summarize these results as 


P 


y,y ’ 


(t) 



_d_ 

dt 


Sino^ (t ) 

y 


(A.114) 


where the approximation of Eq. (A.l4) Is valid whenever Eq. 

(A. 13) Is valid, which requires that the response to the "slow” 
component of turbulence Wg(t) be negligible In comparison with 
the response to the "fast" component Cf(t) z(t), and further- 
more that Eq. (3.^0) be satisfied. 


Let us now consider the ratio ay^/Cy^r. It is well known — 

e.g.. Ref. 11, pp. I 9 O-I 92 of the Wax edition or Ref. 12, p. 48 — 

that /cy! can be expressed In terms of the autocorrelation 
z z 

function (}>y^(T) of ygCt) by 


a 


y 


z 




(0)“ 

^z 

K 

(0) 

^z 



(A. 15 ) 


The right-hand side of Eq . (A. 15) has a simple Interpretation 
that Is analogous to the definition of the Taylor microscale; 
e.g.. Ref. 19j p. 42. Let us approximate <l)y (t) by the para- 
bola <{>y (T)p: 


4> ( T ) - (J) (0 ) 

^y^ p y^ 


1 - 

^ 2 


(A. 16) 
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It follows directly that the parabolic approximation to 4)y^(x) 

given by Eq. (A.l6) becomes zero at t = However, dif- 

ferentiating Eq. (A.l6) twice, we find 




( 0 ) 


Hence, we have 

'*’y_(o) 


(p” (0) 

yz p. 


/2 


= 0.71X, 


(A. 17) 


(A. 18) 


that is, when a parabolic approximation to 4>y^(x) is used, the 

right-hand side of Eq. (A. 15) is slightly less than the time 
delay associated with the zero crossing of (j>y^(x)p. In the tur- 
bulence application, Xq is the analog of the Taylor microscale. 


Equation (A. 15) has another interpretation. It follows 
directly from pp. 192-193 of the Wax edition of Ref. 11 or from 
Eq. (1.65) of Ref. 12 that (cy^/Oyt) is equal to (1 /tt) times the 

expected time between zero crossings of the Gaussian process 
{yz('t))- Now, we generally would expect Xq to be about one-half 
of the nominal "correlation time" of the process. Consequently , 
we may aonalude from both of the above -interpretations that 
(cSu /Oijt) typically is of the order of one-third of the nominal 

"correlation interval" of of the process yz^i^- 


Using the above interpretation of Oy /oyi, we may now inter- 

pret the right-hand side of Eq . (A.l^). Letting the symbol ~ 
denote "is of the order of," we have 


y,y' 


(t) - 


1 

F 


_d_ 

dt 


y 


al ( t ) 


cor 


.2 

I 

y 


ai(t ) 


(A. 19) 


i.e., p ,(t) is of the order of one-sixth of the fractional 

id j ij 2 

change of Oy(t) that occurs in one nominal correlation interval 
of the process yz('l^)- In the present application, the fractional 
change in Oy(t) is approximately equal to the fractional change 
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2 * o 

in Of(t) — see Eq . (A. 7) and recall that Cy is a constant. 

" 2 

Hence, we conclude that when the fractional change in Cf(t) is 
negligible over the nominal correlation Interval of the re- 
sponse process, we have Pv,y*('*^) ~ Finally, we note that 

the requirement that the fractional change in bf(t) over the 
correlation interval of the response process be negligible is 
essentially equivalent to the requirement for the validity of 
the quasi-stationary approximation of Eq. (3-40), as is evident 
from the discussion in Sec. 5 of this report. Thus ^ whenever 
the approximation of Eq. (3,40) is vatidy it is permissible to 
assume that the correlation coefficient between y(t)\o-^ and its 
first derivative is negligibly small. 
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APPENDIX B 

DERIVATION OF CORRECTION TERM TO GAUSSIAN APPROXIMATION 

OF EXCEEDANCE PLOTS 


/ 2 ^ 

To develop an expression for the correction term '^(y|cri) 
In Eq. (^. 36 ), we require, according to Eqs . (4.33) and (4.34), 
an expression for the second derivative with respect to of 
N.(ylai). See Eq . (4.26). According to Eqs. (4.12), (4.14), 
and (4 t 15), N_j_(y|a^) can be expressed as 


N+(y lap 



(B.l) 


where, according to Eqs. (4.l6) and (4.19), we have 






z 


(B.2) 


and 


2 

a. 


y 



(B.3) 


and where the explicit dependence of and on ai has been 

«y «/ 

deleted In the notation of Eq. (B.l) and In the left-hand sides 
of Eqs. (B.2) and (B.3)* 

Forming the first derivative of N^_(y|ap with respect to 
ap we have, using Eq. (B.l) and the notation of Eq. (4.26), 


(y lap 



N^(y lap 


Carrying out 


the differentiation yields 
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- JLl. 
2ot 




[2(op^ da|J 



da? 

da^ -| 

a^ — i - a? 

-JL 

y daj y 

da^ 

(ap^ 


L y 

_ 


da^ 


da: 


da|. 


= I N+(y|ap 


da|. 


o' 


Zl 


y 


da 


2 -1 


da^ 

a? 


(B.4) 


where we have used the definition of Eq. (B.l) 
(B.2) and (B.3)> we have 


But from Eqs. 


. 2 

ZUL 
. 2 
'f 


da; 

h 

da: 


da? 


= o' 


J. = 


da^ 


a: 


(B.5a,b) 


hence, Eq. (B.4) can be expressed as 


N|^\y|ap = I N^(y|ap 


a? 


(B.6) 
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Differentiating Eq. (B.6) with respect to ai yields 


N+^^(ylcrp = I N^(y|ap 


y. 


y‘ 




da 


(<yp^ ddj 


(a|)^ daj 


\ 

^ W 


y. 


-l2 


- 1 . + 


y, 


7 


N+(y |hp 


(B.7) 


Substituting Eqs . (B.5a,b) Into Eq . (B.7), and simplifying the 
resulting expression gives, finally. 


N|^^(y lap 



(B.8) 


VJhen a^|a^ and a^a^ are written for and oi in Eq. (B.8), 

y X y I y ^ 

and the result is combined with Eq. (^.33)> we obtain the 
definition of Q^^^(y|ap given by Eq. (4.34). 
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APPENDIX C 

DERIVATION OF CORRECTION TERM TO GAUSSIAN APPROXIMATION OF 

PROBABILITY DENSITY FUNCTION 


According to Eqs. (^.49) through (4.51) j the second deriva- 
tive of p(ylo|«) with respect to a|. is required to derive an 

expression for the correction term U^^^(y|af) to the Gaussian 

probability density p(ylaf). Suppressing the dependence of Oy 
on- Of, as given by Eq. (B.l), we have from Eq. (4.4) for p(y[cr|'), 

- jlI- 

p(y|ai) = -=^=^ e ^^y . (C.l) 


From the fact that (day/da|>) = Oy [see Eq 


from^Eq. 
to Of is 


(C.l) that the first derivative of 






(B.5a)], it follows 
P(y|cTf) with respect 



(C.2) 


2 • 

Differentiating Eq. (C.2) with respect to and again using 
Eq. (B.5a), we have 
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(We Remind the reader that we have used the notation Oy for 
Oylof.) By combining Eqs. (^t.49) and (C.3), we obtain^ the 

definition of given by Eq. (4.50). 
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APPENDIX D 

DERIVATION OF EXPRESSION FOR COEFFICIENT OF 
EXCESS OF RESPONSE PROCESS 


( 2 ) 

To show that the coefficient of excess Yy of the response 
process y(t) Is given by Eq. (4.5^), we first note that the mean- 
square response E{y^} is obtained from Eq. (4.35a) by 


E{y^> = af, = 


(aMap p(cpda^ 


= + 


p(apdaj 


a 


°y * °f “y 

S Z 


(D.l) 


The fourth moment of the response y(t) may be expressed as 


E{yM = 


E{y‘* lap pCapdCf 


(D.2) 


Substituting Eq. (4.6) into Eq . (D.2), then using the expression 
for a^|a^ given by Eq . (4.35a), and then simplifying, we have 


E{y**} = 3 


(aMop^ p(apdaj 


= 3 


= 3 


s ^ 


+2aJ aj aj. + oj; ap p(apda^ 
s s z z 


= 3(a‘ 




(D.3) 
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f ll-S 

Noting that Is the fourth central moment of y and that the 

J 

mean value of y is zero, we can see from Eqs. (4.52), (D.3), and 
(D.l) that the coefficient of excess of y may be expressed by 


S 2 


^s ^S ^ ^Z 


3 ) 


y 

3(c" )" E{(ap2 _ [E(op]M 
^ z 


y 


(aj )" 

. 3 u'i> 

y 


(D.4) 


which is the expression given by Eq. (4.54), when it is recognized 
that (aMap = 

O' O’ 
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APPENDIX E 
» 

DERIVATION OF EQUATION (5.19) 

To establish the validity of Eq. (5.19)* we shall first 
show that for ergodic processes v'(t), we have 

Var{v”(t) + [v’(t)]2} = E{[v*’(t)]M + Var { [v » ( t ) ] ^ } . (E.l) 

Let us define 

x(t) ^ V'(t) + [v'(t)]" . (E.2) 

Using a well known result, we have 

Var {v ' ' ( t )+[v ’ ( t ) ] ^ } = Var{x} = E{x^} - E^{x} 

= E{[v’ '+(v' )^]M - (E{v* ’ }+E{(v' )2})^ 

= E{(v")M + 2E{v'»(v’)^> + EfCv’)**} 

- E^{v'’} - 2E{v* ’}E{(v’)^> - E^{(v’)M 

= Var{v''} + VarfCv’)^} 

+ 2(E{v'’(v')M - E{v» ’ }E{(v’ )^>) (E.3) 

However, v(t) is a stationary process; hence, E{v’’} = 0. Using 
this result, we may write Eq . (E.3) as 

Var {v ' ’ ( t ) +[v’ ( t ) ] ^ } = E{(v'')^} + Var{(v*)^} + 2(E{v ' ' (v ' ) ^ ) ) • 

(E.4) 

To evaluate the last term in the right-hand side of Eq. (E.4), 
we shall assume that {v'(t)> is an ergodic process; hence, we 
have 
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E{v’ * (v' ) ^> 


lim 


1 

T 


fT/2 

I v"(t)[v'(t)]Mt 
-T/2 


lim 1 

* P ->00 ^ 


Cv’(t)]^ 


T/2 


fT/2 


-2 


V '(t)[v»(t)]2dt 


l-T/2 -t/2 


,(E.5) 


where the second line was obtained by integration by parts. It 
follows from Eq. (E.5) that 


“j v”(t)[v'(t)]2dt = ^ {[v'(T/2)]3 _ [v'(-T/2)]3} . (e.6) 

-T/2 


However, for a stationary process v'Ct), the quantities v'(T/2) 
and v*(-T/2) are bounded; hence, the right-hand side of 
Eq. (E.6) vanishes in the limit T-h». Consequent ly , we have 

E{v* • (v’ = 0. (E.7) 

Therefore, from Eq. (E.4), we have 

Var{v' ' (t)+[v’ (t ) ]^ } = E{[v’’(t)]^} + Var { [v ' ( t ) ] ^ } , (E.8) 

which is the result given by Eq. (E.l). 

Unfortunately, the quantity Var { [v ' ( t ) ] ^ } depends on fourth 
order statistics of the process v'(t). These statistics are 
Impossible to evaluate from the autocorrelation function of 
v(t), unless it is assumed that v'(t) is a (stationary) Gaussian 
process, which necessarily has zero mean value because the 
process {v(tH is, by definition, stationary. If this Gaussian 
assumption is made, it is known (e.g.. Ref. 21, p. 52) that 

Var{[v‘ (t)]2> = 2 e2{[v» ( t)]2) (E.9a) 

= 2[R^*(0)]2 , (E.9b) 

where the second line is a consequence of Eq. (5.15). Further- 
more, from Eq. (5-16), we have 

E{[v» »(t)]2} = R^^^(O). (E.IO) 
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By combining Eqs . (E.8), (E-9)j and (E.IO), we have, finally, 
Var{v'’(t) + [v’(t)]2} = + 2[R^’(0)]", (E.ll) 


which is the result given by Eq . (5.19), and is strictly valid 
only for cases where {v'(t)> is an ergodlc Gaussian process. 
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APPENDIX F 

DERIVATION OF RELATIONSHIP BETWEEN AUTOCORRELATION FUNCTION 
OF A STATIONARY GAUSSIAN PROCESS AND AUTOCORRELATION 
FUNCTION OF THE SQUARE OF ITS LOGARITHM 


Here, we shall derive an expression for the autocorrelation 
function of An z^(t); l.e.. 


r'(t.) = E{An z^(t) An z^(t+x)} (F.l) 

in terms of the autocorrelation coefficient of the stationary 
Gaussian process z-^(t) , where 

E{Zj^} = 0 , E{zJ} = 1 ; (P.2a,b) 


l.e.. 


P(t) = E{Zj^(t) Zj^(t + T)} 


(F.3) 


We shall carry out the derivation using the form of Price's 
Theorem given by Eqs . (20) and (21) of Ref. 22. 

Let Zh(t) be a stationary Gaussian process satisfying Eqs. 
(P.2a,b), and let f[zh] be an arbitrary zero-memory (generally 
nonlinear) transformation of z^. The form of the theorem that 
we shall use states that 


1^ = E{f[z^(t)l f'[Zj^(t+T)]} , (F.il) 

where the primes denote the derivative of fCz^] with respect to 
zyi, and where these derivatives are evaluated at times t and 
t + T, as Indicated. In the left-hand side, p denotes the cor- 
relation coefficient defined by Eq. (F.3) and R denotes 

R(t) = E{f[Zj^(t)] f[zj^(t+x)]} . (F.5) 

Thus, comparing Eqs. (F.l) and (P.5)j It is evident that we are 
Interested in the case 
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9 


(F.6) 


= An 



hence, we have 


f » [z 


h 



(P.7) 


Substitution of Eq. (P.7) Into Eq. (P.4) yields, when the ex- 
pression for the joint Gaussian density Is written out and where 
we substitute Zj = Zj^(t) and = Zj^(t+x), 


3R 

3p 


(zf+Z2“2pZjZ2) 


2tt/1-p ' 




2(1-pM 


dz , dz, 
1 2 


_00 — .00 


(P.8) 


Let us define 


yi = 


/2(l-p2) 


Y2 = 


/2(l-p2) 


Then, Eq . (P.8) reduces to 


3R ^ 

3p 


2 

7r/l-p^ 


fOO 


j 

«00 



1 

yjy2 


-(yi+y|-2pyiy2 ) 

e 


dyidyz 


(F.9a,b) 


(P.IO) 


Using the Integral given on the bottom of p. 207 and the top of 
p. 208 of the Wax edition of Ref. 11, we have, formally, for n 
and m both equal to -1 and p equal to -cos(j>. 


3R 

9p 


r(i-|) r(i-|) 


7r/l-p‘ 


( slncj) ) 


- 1 


cosc})^?^ (1,1;|-; cos^cj)) 


= ±4p2Fi(l,l;|; p") 


(F.ll) 


where we have substituted r(l/2) = /ir and p = -cos(f) (hence, 
sine}) = ±/l-p^) and where ^PiC***) is the hypergeometrlc function. 
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However, using Eqs . (A. 1.35) and (A. 1.39a) of Ref. 23 
1076 and 1077 , we have 

= (I_p2)-H arcslnp ^ 

Thus, combining Eqs. (F.ll) and (P.12) yields 


on pp . 


(F.12) 

I 


M 

3p 


/i-p^ 


arcsinp 


(F.13) 


To determine R, we may Integrate Eq. (P.13) with respect to p: 

■p(t) 


R(t) = ±4 


arcsing 


/1-g- 


dg + c , 


(F.14) 


where c is the constant of integration and g is a "dummy variable." 
Now, when p(t) = 0, Zh(t) and Zh(t+r) are uncorrelated according 
to Eq. (F.3). In this case, it follows from Eq. (P.l) that 


R = {E[5-n = 4{E[£n z^]}^ , for p = 0 , 


h- 


(P.I5) 


which, according to Eq . (F.l4), is our constant of integration. 
Thus, we have from Eqs. (F.l4) and (P.15), 




fp(T) 


R(t) = 4’E^[£n 




arcsing 


dg 


(F.I6) 


However, 


dg 


arcsing = 


/1-g' 


(P.17) 
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therefore. If we let F(5) = arcsln^, the integral in Eq. (F.l6) 
is of the form 

■p(t) p(t) /-p(t) 

F(?) F’(?)dC = [F(C)]^ F(?) F'(?)d5 , (F.18) 

J 

0 0 0 

where the prime denotes differentiation. Thus, from Eq. (F.l8), 
we have 

•p(t) 

F(?) F»(C)d^ = I {[F(p)]2 - [F(0)]M . (F.19) 

0 

Using F(C) = arcsinC and noting that arcsinO = 0, we have, by 
combining Eqs. (F.l6) and (F.19), 

R(t) = 4{E^[Jln Zj^] + ^ arcsin^ p(t)} . (F.20) 

We may now determine the correct sign in Eq . (F.20). When p = 0, 
arcsin p = 0. On the other hand, when p = 1, arcsin p = 7 t/ 2 . 
However, p = 1 occurs when t = 0, and for this value of t, R(t) 
must achieve a maximum. This is possible only with the plus 
sign in Eq. (P.20). Consequently, the correct sign in Eq . (F.20) 
yields 

R(t) = {E[Jln + 2 arcsin^ p(x) . (P.21) 

The result of Eq . (P.21) can be checked as follows. At 
T = 0, we have p(x) = 1; hence, arcsin p = ir/2 at this point. 
Consequently, Eq . (F.21) gives, for this value of t = 0, 

R(0) = {E[£n z^J}^ + ^ . (F.22) 

But, from Eq. (F.l), we have 

R(0) = E{[Jln z^]2} . (F.23) 

By combining Eqs. (F.22) and (F.23), we can see that Eq . (F.21) 
yields the result 
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E{[Jln - {ECiln ^ 


(F.24) 


that is, Eq. (P.21) predicts that the variance of An Is equal 
to it^/2. Since z^ Is, by assumption, Gaussian with zero mean 
and unit variance, the result of Eq. (P.24) can be checked 
directly. Por the expected value of An z^, we have, using the 
fact that An z^ Is an even function of z^. 


E{An z^} 


,00 


1 

/2it 


An z^ e 
h 


~CD 



4 

/2'ir 


•oo 

An z, e 



dz, 

h 


= - (C+An 2) 


(P.25) 


where C Is Euler's constant. 


C = 0.577215... , 


(P.26) 


and where the last line In Eq. (F.25) was obtained using Formula 
(4.333) on p. 574 of Ref. 24. Por the expected value of 
[An z^]^, we have, again using the fact that An z^ is an even 
function of 


E{[An zJ]M 


1 

/2tt 


,00 


(An z^)^ 

_ 0 O 


e 



dZh 


8 

/2it 


• 00 

(An z^)^ e 



0 


(P.27) 


To put Eq. 
substitute 


(P.27) into the form of a tabulated integral, we 

= Zh/2. With this substitution, Eq. (F.27) becomes 
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L 



_ 8 _ 



0 


[(iln /2)" 


+ 2(Hn /2)(lnK) + (AnO^] 



= Ii + I2 + 


(P.28) 


where Ii, I 2 , and I3 are the three terms that result from Inte- 
grating each term within the brackets in Eq . (P.28) separately 

and then multiplying each by the common factor 8//F. Using, from 
Ref. 24, pp. 307 and 574, Formulas 3.321.3, 4.333, and 4.335.2 
to evaluate Ij , and I3, respectively, we find 

Ij = 4(iln /2)^ = (£n 2)^ , (P.29) 

12 = - 4(iin /2) (C + iln 4) = - 2 ( An 2) (C + ^n 4 ) , (P.30) 

13 = (C + 2Zn 2)^ + ^ = (C + £n 4) + ^ . (F.3D 


Combining Eqs. (F.28) through (F. 3 I), we have 

E{[£n z2]2} = 1 ^+ 12+13 

= (£n 2)^ - 2(£n 2) (C + £n 4) + (C + £ri 4)^+ 

=[(£n2)-(C+£n4)]2+^ 

= [(£n 2) - C - 2(£n 2)]^ + ^ 

= (C + £n 2)^ + ^ . (P. 32 ) 
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I 


Finally, combining Eqs. (P.32) and (P.25) we obtain the desired 
result 


E{[iln - {E[An z*]} 



(P.33) 


which is in perfect agreement with the result of Eq. (F.24). 

This completes the check of Eq. (P.24). 

To obtain the final result for R(t), we combine Eqs. (P.21) 
and (F.25); this gives 


R(t) = (C + £n 2)^ + 2 arcsin^ p(t) , (P.34) 

which is the final result. Insofar as we are aware, the result 
of Eq. (F.34) is new. We remind the reader that R(t) and p(t) 
are defined by Eqs. (P.l) and (F.3), and that {zh(t)} is a 
stationary Gaussian process with zero mean and unit variance as 
indicated by Eqs. (F.2a,b). 
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APPENDIX G 

METHOD FOR ESTIMATION OF E{ap 

In Eq. (6.13), a high-pass filtered version 

w^(t) = A' a^(t'> (G.l) 

of the turbulence record w(t) was considered, where A’ z^(t) Is 
the filtered version of the original component z(t) of Eq . (2.3) 
which satisfies 

E{z^} = 1 . (G.2) 


Let Hh(f) denote the high-pass filter complex frequency response 
function used in obtaining wj^(t) from w(t). Then, using the 
approximation of Eq. (6.l8), we may express the mean square 
value of W}^(t) as 


E{»p 


E{o|) 






(0.3) 


where 'I>z(f) is the power spectral density of the component z(t), 
which satisfies, according to Eq. (G.2), 

,oo 

$ (f)df = 1 . (G.4) 

. 

«.oo 


Solving Eq. (G.3) for E{of} yields 

E{wJ} 

E{ap = . (G.5) 

I |H^(f)|2 <D^(f)df 
— 00 


According to the discussion in Sec. 6.1, $z(f) may be assumed 
to have the appropriate von Karman form; thus, because $ 7 ;(f) 
is constrained by Eq. (G.4), it is defined by its integral scale. 
A method for estimating the integral scale was suggested in 
Sec. 6.1. Furthermore, lHj^(f)|* is a known function and 
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E{w^} can be measured from the filtered waveform w^Ct). There- 
fore, Eq. (G.5) can be used to estimate E{a|*}. 

Let us now turn to Justification of Eq. (6.92). Taking 
the nth moment of both sides of Eqs. (2.3) and (G.l), we obtain 

E{w^> = E{a5> E{z’^} (G.6) 


and 

E{w[J} = (A')*^ E{aJ} E{z[J} , (G.7) 

from which we obtain 


E{w5> E{z"} 

— = . (G.8) 

E{w^} (A’)"^ Eizp 

However, since z(t) and z^(t) are both normally distributed with 
zero mean values, we have 

E{z’^} = E{zJJ} , (G.9) 


where B is a constant. Combining Eqs. (G.8) and (G.9), we have 


E 




(G.IO) 


However, from Eq. (2.3), we can write 

E{wp 

E{ap = , (G.ll) 


where the denominator in Eq . (G.ll) is unity. Combining Eqs. 
(G.5) and (G.ll) yields 
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1 

E{wp = E{wM „ 


which Is of the form of Eq. (G.IO) for n = 2. Equation (6 
follows directly from Eqs. (G.IO) and (G.12), If we define 


K 





-th 


|Hh(f)|^ 


(G.12) 


90) 


(G.13) 



APPENDIX H 

DERIVATION OF GENERAL EXPRESSION FOR MOMENTS OF w.(t) 
IN TERMS OF MOMENTS OF w(t) and w^(t) ^ 


Let w and be independent random variables and denote 

their sum By w: 

w=w+w„ (H.l) 

s I 

Let p , p and p denote the probability density .functions of 

Vf V/ V? rt 

S 1 V , 

w, w , and w„. It is well known — e.g.. Ref. 9 — that p^ is the 
convSlution of p^ and p^ ; i.'e., 

f 


P (w) 


P (w-C)p (?)d^ 
s f 


(H. 2a) 



(5)p (w-?)dC. 

Wf 


(H.2b) 


Taking the nth moment of both sides of Eq. (H.2b), we have 


M 


(n) A 


w 


w^p^(w)dw 



.cx> _00 


(?)p (w-C)dCdw 

Wf 


fPw 


.00 

w^p (w-C)dw 

W 

-OO 


I 

00 


d? 




I 

..oo 

J (5+u) p^ (u)du 

oo ^ 


d? 


where we have Introduced the change of variable u = w 
Expanding (5+u)^, we have 
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(n) 

w 



(?) 





k n-k) 
u \ 


>p (u)du d? 

) Wf 



(?)dC 



*00 


(u)du , 


(H.3) 


where we have used the binomial expansion 


(?+u)’^ 



k n-k 
u , 


where 


(H.4) 



n! 

k! (n-k) ! 


(H.5) 


are the binomial coefficients. 

The last line In Eq. (H.3) is a relationship among the 
moments Of p p 

s '^f 


" k=0^*' “s "f 


k=0 ^ ”f 


(H.6) 


Since the area under a probability density Is unity, by defini- 
tion, we always have 


W W^ 


= = 1 
^S 


(H.7) 


Using this fact, we may write out the relations of Eq. (H.6) 
for n = 1, 2, 3> ^ as 
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M 




w 


w. 


w. 


M 


w w^ w« w_ w_ 




p^(3) ^ j^(3) + 3 iv[( 2) (1) gj^Xl)j^^(2) _^. jyj(3) 

W W^ Wg Wj. Wg 


w w~ w~ w„ w~ w„ w- w„ w„ 

I rs IS IS s 


The above relationships may be solved successively for M 


( 1 ) 


f 2) f 

w ’ w ’** yield the set of relationships given by 
s s 

Eqs. (6.89) in the main text. The general form of these 
relationships Is easily seen to be given by Eq. (6.87) 
i . e . , 


"w"’ = ”w"’ - ‘'’''w*'’- " ' 

'^s ^ k=0 ^ '^f ^s 


(H.8) 


(H.9) 
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